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Section  1.  - MODELLING 


1 e 0 Introduction 


I 


This  section  describes  the  results  of  the  first  part  of  the  NAVSTAR/ 
GPS  Navigation  Analysis  and  Algorithm  Development  Study.  The  purpose 
of  this  first  sUption  is  to  define  the  models  for  use  throughout  the  entirety 
of  the  study.  *fhe  ultimate  objective  of  the  study  is  to  determine  a set  of 
algorithms  which  can  be  used  in  a NAVSTAR/GPS  user  navigation  system. 
The  criteria  for  the  acceptability  of  the  algorithm  will  be  the  accuracy  of 
the  position  and  velocity  determination.  Since  there  is  no  way  of  getting 
real  data  for  algorithm  verification  at  this  time,  it  is  imperative  to  estab- 
lish system  and  error  models.  The  results  of  analysis  of  proposed  algo- 
rithms is  then  relative  to  the  models  used. 

An  additional  objective  of  the  entire  study  is  that  it  will  be  a design 
aid.  Thus  it  is  intended  that  through  analysis  of  different  receiver  config- 
urations, an  analytic  basis  for  certain  design  decisions  can  be  established. 

y 

The  final  computer  program  for  analysis  will  not  be  an  interactive  "com- 
o/ 

puter  aided  design^ tool)  however,  it  will  provide  for  inclusion  or  exclusion 
of  certain  receiver  options  and  variability  of  certain  parameters  in  order 
to  make  the  analysis  useful  for  making  design  decisions.  te*- 

The  modelling  effort  for  this  task  is  divided  into  two\«>arts,  Sys- 
tem Model  and  Error  Model.  The  System  Model  contains  the  satellite 
model,  user  model,  control  system  model,  and  data  stream  model.  In 
addition,  scenarios  are  defined  for  evaluation  of  the  proposed  numerical 
algorithms.  The  Error  Model  defines  the  contributors  to  the  navigation 
error.  This  includes  errors  in  the  system  components,  such  as  receiv- 
er and  transmitter  errors,  and  error  sources  exogenous  to  the  System 
Model,  such  as  atmospheric  effects.  The  distinction  between  System 
Model  and  Error  Model  is  somewhat  arbitrary  in  some  instances.  How- 
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ever,  since  the  models  are  not  intended  to  be  used  separately,  this  will 
not  cause  any  difficulty.  The  emphasis  in  this  report  is  the  form  of  the 
models.  Specific  parameter  values  may  be  changed  in  the  course  of  the 
study;  however,  the  form  of  the  models  should  remain  constant. 

The  parts  of  this  report  dealing  with  the  system  model  have  been 
taken  from  various  references  (1,  3,  4,  6,  9,  15)  and  the  applicable  parts 
selected  for  inclusion.  A critical  look  has  been  taken  at  each  model,  but 
no  new  models  are  presented.  The  error  models  are  taken  from  various 
references  (7,  10,  11,  12,  14)  also  along  with  standard  linear 'models  for 
certain  error  types.  Each  error  source  was  examined  and  the  appropriate 
model  chosen  according  to  accepted  modelling  procedures. 

1.  1 System  Model 

The  System  Model  must  be  defined  as  completely  as  possible  at 
the  outset  of  the  study  since  the  model  serves  as  a set  of  ground  rules. 

The  only  major  component  of  the  system  which  will  not  be  modelled  is 
the  computational  unit.  This  will  have  an  effect  on  the  algorithms;  how- 
ever, this  effect  is  clearly  beyond  the  scope  of  this  study.  The  only  config- 
uration unknowns  in  the  subsystems  which  are  modelled  are  in  the  area  of 
the  receiver  model.  Modelling  of  various  types  of  user  receivers  will 
aid  in  the  development  of  the  design  specifications  for  the  receiver.  The 
System  Model  is  divided  below  into  several  segments  solely  for  the  pur- 
pose of  exposition. 

1.  1.  1 Ephemeris  Model 

^e  ephemeris  model  is  concerned  with  the  satellite  constella- 
tion,  the  mathematical  model  for  satellite  motion,  and  the  information 
about  the  satellites  which  the  control  system  segment  supplies  to  the 
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1.  1.  1.  1 Satellite  Constellation 


The  Phase  III  satellite  constellation  consists  of  twenty-four  satel- 
lites. The  constellation  will  have  three  planes  of  satellites  in  approxi- 
mately circular  twelve-hour  orbits.  The  orbit  planes  are  inclined  at  ap- 
proximately 63°  and  spaced  so  that  the  ascending  nodes  of  the  orbit  planes 
are  120°  apart.  Each  orbit  plane  has  eight  satellites  equally  spaced  in 
the  orbit.  It  is  intended  that  between  6 and  11  satellites  will  be  visible 
from  any  point  on  the  Earth  at  all  times.  On  the  average  there  will  be 
eight  or  nine  satellites  in  view.  A complete  definition  of  the  Phase  III 
satellite  constellation  can  be  found  in  reference  1. 

1.1. 1.2  Satellite  Motion 

The  satellite  motion  generation  for  the  System  Model  will  use  a 
simple  two-body  orbit  for  each  satellite.  The  positions  of  each  satellite 
with  respect  to  the  Earth  could  then  be  computed  as  a function  of  time 
from  the  six  orbital  elements  (2),  The  model  used  by  the  system  for  pre- 
diction of  satellite  position  uses  fourteen  elements  to  achieve  the  desired 
navigational  accuracy  (3).  Linear  perturbations  to  this  fourteen  element 
model  can  be  expressed  by  the  orbital  element  model.  Therefore  to 
conserve  computation  in  the  error  analysis,  this  simpler  model  can  be  used. 

The  nominal  GPS  orbit  configuration  calls  for  circular  orbits.  For 
analysis  purposes,  this  simplifies  the  satellite  position  and  velocity 
calculations  since  only  four  parameters  must  be  specified;  the  radius 
of  the  orbit,  orbit  inclination,  longitude  of  ascending  node,  and  time 
of  passage  of  the  ascending  node.  In  fact,  all  of  the  satellite  locations 
and  velocities  can  be  computed  by  orthogonal  transformations  of  a 
single  satellite  position  and  velocity 

Xi  = TiXl  1 = 2»  3»  •••»  24 
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•1 


where  X.  is  the  state  vector  for  the  i satellite  and  T,  is  an  orthogonal 
1 1 

(length  preserving)  transformation. 


1.  1.  1.  3 Data  Link  Information 

The  control  system  segment  of  GPS  as  a part  of  its  function  sup- 
plies orbit  information  to  the  user.  There  are  two  types  of  orbit  data 
provided.  One  type  is  the  very  accurate  set  of  fourteen  parameters  which 
are  updated  every  day,  the  other  type  is  the  Almanac  data. 

The  very  precise  data  consists  of  fourteen  parameters  which  are 
to  be  used  in  conjunction  with  the  nominal  orbit  parameters.  The  satel- 
lite is  updated  once  a day  with  twenty-four  sets  of  this  data  Each  data 
set  is  optimized  to  fit  the  predicted  orbit  accurately  over  one  hour.  The 
data  sets  should  not  be  considered  as  true  orbit  parameters,  but  rather 
as  coefficients  of  a numerical  fit  to  the  predicted  orbit.  The  parameters 
take  the  form  of  orbital  elements  and  additional  information  to  account 
for  some  of  the  error  introduced  by  effects  not  considered  in  the  standard 
two-body  orbit  (viz.  pole  wobble,  gravitational  anomalies,  solar  pres- 
sure, etc. ).  These  corrections  are  valid  only  for  the  instant  at  which 
they  are  computed  and  tend  to  degrade  with  time.  Consequently,  the  con- 
trol system  segment  tracks  the  satellites  and  updates  these  parameters 
on  twenty-four  hour  intervals  (3).  The  most  current  information  available 
(4)  lists  the  parameters  as: 


M - mean  anomaly  at  reference  time 
o 

An  - mean  motion  deviation 


- eccentricity 

- square  root  of  semi-major  axis 
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0 - right  ascension  at  reference 
o 

1 - inclination  at  reference 

o 

u)  . argument  of  perigee 

e 

0 - rate  of  right  ascension 

C - amplitude  of  the  cosine  harmonic  correction  term 
to  the  argument  of  latitude 

C - amplitude  of  the  sine  harmonic  correction  term 
to  the  argument  of  latitude 

C - amplitude  of  the  cosine  harmonic  correction  term 
to  the  orbit  radius 

C - amplitude  of  the  sine  harmonic  correction  term 
to  the  orbit  radius 

C.  - amplitude  of  the  cosine  harmonic  correction  term 
to  the  angle  of  inclination 

C.  - amplitude  of  the  sine  harmonic  correction  term 
to  the  angle  of  inclination 

In  addition,  the  following  parameters  are  provided. 


t - reference  time  ephemeris 
oe 

AODE  - age  of  data  (ephemeris) 


The  algorithm  to  convert  this 
as  follows: 

fi  = 3.986008  X 1014  --te™- 

sec 

O = 7.292115147  x 10'5  — 

e sec 

x 

A = (Ya)2 


data  to  satellite  position  in  ECI  is 

WGS  72  value  of  the  Earth's  uni- 
versal gravitational  parameter 

WGS  72  value  of  the  Earth's  rota- 
tion rate 

semi-major  axis 


j 
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n = 

° V A3 


computed  mean  motion 


t , = t - t * 
k oe 


time  from  epoch 


n = n + An 
o 

M = M + nt, 
k o k 

M = E - e sin  E, 
k k k 

cos  v = (cos  E - e)/(l-e  cos  E.  ) 


sin  vw  = V1  - e sin  E /(1-e  cos  E,  ) 


corrected  mean  motion 

mean  anomaly 

Kepler's  equation  for 
eccentric  anomaly 

true  anomaly 


•k  = vk + " 


6uk 

= C 

us 

sin  26. 
k 

+ C 

uc 

cos  26^ 

6rk 

= C 

rc 

cos  26, 
k 

+ C 

ri 

, *ln  2*i 

tl  ■■ 

= c. 

cos  26. 

+ c. 

sin  26, 

k 

1C 

k 

is 

k 

Uk  = 

*k  + 

fiuk 

rk  = 

A<1  . 

e cos  E 

k»  + # 

rk 

ik  = l,'  6ik 


argument  of  latitude 

> 

argument  of  latitude 
correction 

second 

radius  correction  “harmonic 

perturbations 

correction  to  inclina- 
tion j 

corrected  argument 
of  latitude 

corrected  radius 
corrected  inclination 


x.  = r.  cos  u, 
1 k k 


yk  = rk  8m  “k 


0.  = 0 + (0  - 0 )t.  - o t 
k o e k e oe 


position  in  orbital 

plane  / 

/ 

corrected  longitude  of 
ascending  node 


t is  GPS  system  time  at  time  of  transmission,  i.e.,  GPS  time  corrected 
for  transit  time  (range /speed  of  light) 
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r 


*k  * co’  °k  • yk  co*  Sc  ,in  °k 

yk  - .in  Ok  + y;  co.  ik  co. 

*k  * yk  ,ln  S. 

Tentative  scaling  and  resolution  information  for  each  parameter 
are  contained  in  reference  4. 

The  data  stream  from  each  satellite  includes  Almanac  data  for 
all  of  the  satellites.  This  data  consists  of  orbital  elements,  satellite 
ID  and  health,  and  time  parameters  (4).  The  Almanac  parameters  may 
be  used  for  alert  calculations  to  determine  which  satellites  are  "in  view", 
i.  e.  , which  may  be  received.  It  is  not  intended  that  Almanac  data  be 
used  for  precise  navigation,  but  it  will  be  useful  for  acquisition.  Al- 
manac data  is  updated  every  six  days  (4). 

1.1.2  Time  Model 

The  time  model  is  concerned  with  the  satellite  clock  and  deter- 
ministic delays  in  the  downlink.  The  system  model  of  time  (3,4)  has 
terms  in  it  so  that  the  user  can  correct  the  measured  time  of  arrival 
of  the  satellite  signal  for  the  deterministic  part  of  clock  off-sets,  biases, 
satellite  electronic  delays,  relativistic  effects,  and  atmospheric  delays. 

1. 1.2.1  Clock  Model 

The  satellite  oscillator  frequency  is  set  at  a nominal  value  which 
includes  a frequency  offset  to  account  for  the  relativistic  effects  in  the 
nominal  orbit.  The  control  system  segment  of  GPS  then  monitors  each 
satellite  vehicle  clock  to  calibrate  deterministic  errors  in  the  present 
value  such  as  offset,  frequency  bias,  and  relativistic  effects  due  to  off- 
nominal  orbits.  The  control  system  then  computes  three  parameters 
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with  which  the  user  may  correct  the  system  time.  These  are  sent  as 
part  of  the  data  from  the  satellite  to  the  user.  The  correction  param- 
eters are  the  coefficients  of  a polynomial  correction  which  the  user  can 
apply.  This  correction  has  the  form: 


At  = art+ 


a,  (t 
1 s 


toc) 


a2(ts  ' 


t ) 
oc 


so  that 

t = t + 

8 

where 

a^,  a^,  a^  are  the  correction  parameters 

t is  the  satellite  vehicle  time 
s 

t is  the  satellite  vehicle  clock  epoch  time 
oc  \ 

t is  the  corrected  system  time 

The  coefficients  a^,  a^,  and  a^  are  computed  to  fit  the  predicted 
clock  behavior  over  small  time  intervals.  The  control  system  sends  new 
sets  of  these  parameters  to  the  satellite  every  day. 

For  t_toc  < min* » t^ie  approximations  provide  errors  less  than 
. 5 nanosecond.  After  45  minutes,  these  errors  of  approximation  degrade 
as: 

T-toc  error 


1 hr 

1 ns 

2 hr 

8 ns 

3 hr 

26  ns 

4 hr 

60  ns 

8 
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This  polynomial  form  of  the  correction  does  not  provide  for  a 
graceful  degradation  of  the  relativistic  errors.  A more  graceful  degra- 
dation of  error  can  be  achieved  by  using  the  following  correction  (4): 


At  = (a 


a0r,+ 


(al  ‘ 


alr,(t. 


•oc*' 


+ (a. 


a ,J(t  - t ) 
2r  s oc 


+ At 


where 


K ■ (■ 


4.  443  X 1 0 


-10  sec  \ 


eV^.inE  (t) 


aQr  = -4.443  X IQ’10  e VX  sin  E (t  ) 


air  = _4*443  X 10"10  tfSm  eVS  n COB  E (toc)/[1  - e cos  E 

a2r  = 2,2215  x 10  vffKBSr  eV^n2  sin  E (tog)/!1  - e cos  E (t^)]2 


The  rationale  for  this  type  of  correction  can  be  found  in  refer- 
ence 14  and  in  references  18  and  52  of  reference  14. 

1«1.2.2  Atmospheric  Model 

The  atmospheric  model  is  constructed  to  predict  the  range  error 
introduced  by  propagation  error  due  to  the  atmosphere,  specifically  the 
ionospheric  and  tropospheric  delay.  Of  the  two,  the  tropospheric  error 
has  less  variation  and  the  range  error,  AR,  can  be  approximated  to  bet- 
ter than  1 ft  at  elevations  greater  than  10°  and  better  than  5 ft  at  eleva- 
tions between  5°  and  10°  by  the  simple  formula  (3) 

AR  = K esc  E 
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or  (12) 


AR 


K 


1 

sin  E + . 026 


where 


E,  the  elevation  angle  from  the  user  to  the  satellite  is  greater 
than  5°  and  from  (3)  K is  a constant  (•■  4 ft). 


The  user  will  use  one  of  these  models  to  compensate  for  tropospheric 
range  errors. 

The  ionospheric  model  must  account  for  seasonal  and  diurnal 
variations  and  latitude  dependencies  (3).  There  are  two  correction 
schemes  available  to  the  user.  The  first  method  uses  eight  parameters 
which  are  sent  in  the  data  stream.  It  should  be  noted  that  the  following 
model  is  new  and  documentation  as  to  its  validity  is  not  yet  available. 

The  ionospheric  correction  time  is  calculated  by 

IONO 


F* 


IONO 


5.  x io*9  + ( £ a <pn 
\n=0  n m 


F*  5.  X 10 


-9 


where 


x = 


(t  - 50400) 
Sjl  (Pn 

n=0n  m 


10 


. 57 
(sec) 
. 57 


M. 


and 


F = 1.  + 16.  [.  53  - E]" 


a and  fi  ; n = 0,  1,2,  and  3;  are  the  satellite  transmitted  data  words, 
n n 


Other  equations  that  must  be  solved  are 


t = 4.  32  X 10  Xj  + GPS  time  (sec);  t > 86400  use  t = t - 86400 


<p  = <p.  + 0.064  cos  (X.  - 1.617)  (semi-circles) 

mi  i 


X.  = X 4 ^ (semi-circles) 

1 u cos. 


*i  = 


(Pu  + ^cos  A (semi-circles),  <p  < .416  (semi-circles) 


<P 


u 


, <Pu  > .416  (semi-circles) 


0.0137  _ . „ 

♦ = £~~'q  j j - 0.022  (semi-circles) 


The  terms  used  in  computation  of  ionospheric  delay  are  as  follows: 

* Satellite  Transmitted  Terms 

0tn  - the  coefficients  of  a cubic  equation  representing  the 

amplitude  of  the  vertical  delay  (4  coefficients  - 8 bits 
each) 

Pn  - the  coefficients  of  a cubic  equation  representing  the 
normalised  period  of  the  model.  The  true  period 
has  been  divided  by  2fT.  (4  coefficients  - 8 bits  each) 
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• Receiver  Generated  Terms 


E - elevation  angle  between  the  user  and  satellite  (semi- 
circles) \ 

A - azimuth  angle  between  the  user  and  satellite,  mea- 
sured clockwise  positive  from  the  true  North  (semi- 
circles) 

<p  - user  geodetic  latitude  (semi-circles)  WGS-72 

X - user  geodetic  longitude  (semi-circles)  WGS-72 
u 

GPS  time  - receiver  computed  system  time 
Computed  Terms 

F - obliquity  factor  (dimensionless) 
t - local  time  (sec) 

O - geomagnetic  latitude  of  the  Earth  projection  of  the 
m 

ionospheric  intersection  point  (mean  ionospheric 
height  assumed  350  km)  (semi-circles) 

X.  - geodetic  longitude  of  the  Earth  projection  of  the  iono- 
spheric intersection  point  (semi-circles) 

O.  - geodetic  latitude  of  the  Earth  projection  of  the  iono- 
spheric intersection  point  (semi-circles) 

^ - Earth's  central  angle  between  user  position  and  Earth 
projection  of  ionospheric  intersection  point,  (semi- 
circles) 

The  values  of  0r0,  Otj,  Of 2,  and  oty  and,  fiQ,  fly  and  03  are  transmitted 
in  Data  Block  I with  8 bits  /coefficient  or  64  bits  total. 

For  a dual  frequency  receiver,  an  alternate  correction  scheme  may 
be  used  to  calculate  the  user  range  error.  The  range  error  is  calculated 
according  to  (3) 
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» 


> where 

Rj  is  pseudo- range  at  frequency  f^  (denoted  L^) 

R2  is  pseudo-range  at  frequency  f2  (denoted  L 

^Rj  is  the  range  error  of  R^ 

The  dual  frequency  ionospheric  delay  correction  scheme  does  not 
require  any  external  data  source.  It  does  require  reception  of  two  dis- 
tinct frequencies,  Lj  and  L^,  from  the  satellite.  This  correction  scheme 
removes  more  of  the  uncertainty  due  to  ionospheric  delays. 

1.1.3  Receiver  Model 

> 

One  of  the  purposes  of  this  entire  study  is  to  aid  in  the  design  of 
the  receiver  by  performing  trade-off  studies.  For  this  reason,  the  base- 
line receiver  model  has  not  been  established.  This  section  will  thus  present 
the  alternatives  for  the  trade-off  studies  and  point  out  other  areas  of  the 
receiver  model  which  impact  the  navigation  problem.  The  trade-off  areas 
are  1)  single  channel  or  multiple  channel  receiver,  2)  range  and  range 
rate  or  range  only  data,  and  3)  whether  to  demodulate  the  incoming  infor- 
mation stream.  Other  areas  of  importance  to  the  navigation  problem  in- 
clude the  data  interface  between  the  receiver  and  the  processor,  acquisition 
time,  and  the  frequency  of  measurements.  The  accuracy  of  the  receiver 
will  be  discussed  in  the  error  model  section. 

At  the  present  time,  there  are  three  candidate  receiver  configura- 
tions. They  are  1)  dual-frequency  multi-channel,  2)  single-frequency  sin- 
gle channel  for  P code  and  C/A,  and  3)  single -frequency  single  channel 
C/A  only. 
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1.  1.  3.  1 Receiver  Channels 


The  trade-off  in  the  number  of  channels  is,  of  course,  cost  vs. 
navigational  accuracy.  In  addition,  the  time  to  first  fix  and  through  put 
must  be  considered.  The  following  describe  the  processing  considera- 
tions of  single-channel  and  four-channel  receivers. 

1.  1.  3.  1.  1 Single  Channel  Receiver 

In  order  for  a single -channel  receiver  to  supply  sufficient  data 
for  accurate  navigation,  it  must  switch  between  satellites  for  each  new 
piece  of  data.  This  switching  involves  the  acquisition  time  for  each  new 
satellite.  This  is  particularly  important  in  a single  fix  case  where  four 
different  satellites  must  be  received  before  any  position  computation  can 
be  made.  If  recursive  filtering  of  data  is  to  be  done,  then  each  satellite 
measurement  can  be  incorporated  in  a sequential  fashion  as  it  is  received. 
The  acquisition  time  from  one  satellite  to  the  next  thus  becomes  an  important 
parameter  of  the  navigation  accuracy  in  the  single- channel  receiver. 

1.1.  3. 1.2  Four-Channel  Receiver 

This  is  a special  case  of  a multiple-channel  receiver.  A four- 
channel  receiver  is  being  considered  since  this  is  the  minimum  needed 
to  get  a single  fix  without  having  to  switch  satellites.  With  four  chan- 
nels of  data  reception,  acquisition  time  does  not  effect  the  accuracy  of 
the  navigation  algorithms.  Even  in  a recursive  filtering  mode,  some 
data  is  available  while  other  satellites  are  being  acquired.  With  a four- 
channel  receiver,  either  sequential  or  batch  processing  of  data  may  be 
utilized  and  should  be  analyzed. 


1. 1.3.2  Measurable  a 


\ 


The  trade-off  to  be  made  in  this  area  is  whether  or  not  to  make 
Doppler  measurements  in  addition  to  time* of- arrival  measurements. 
Making  Doppler  measurements  adds  to  the  complexity  of  both  the  hard- 
ware and  the  navigational  software.  The  additional  cost  will  be  traded- 
off  against  the  navigational  accuracy  obtainable.  The  Doppler  measure- 
ments will  be  particularly  useful  in  reducing  platform  position  and 
velocity  uncertainties  when  platform  accelerations  are  taking  place. 

1.1.3.  3 Data  Demodulation 

This  trade-off  concerns  whether  the  ephemeris  and  clock  data 
will  be  extracted  from  the  satellite  signal.  If  the  data  is  not  demodu- 
lated, the  nominal  values  must  be  used  for  these  parameters.  It  is 
possible  that  some  method  for  updating  the  nominal  values  on  some  in- 
frequent basis  can  be  arranged.  Both  the  hardware  and  the  processing 
software  can  be  simplified  if  it  is  not  necessary  to  demodulate  data.  The 
trade-off  here  is  a rather  large  degradation  of  navigational  accuracy. 
(Note  - It  may  be  possible  to  reduce  the  effect  of  not  demodulating  data 
through  increased  computational  complexity.  The  study  of  this  case  is 
beyond  the  scope  of  this  study.  ) 

1.1.  3.  4 Receiver  Frequencies 

The  satellites  transmit  on  two  separate  frequencies  denoted 
and  Lg.  All  of  the  information  required  for  navigation  is  contained  on 
each  frequency.  Reception  on  two  frequencies,  in  addition  to  two  sets 
of  data,  allows  a more  accurate  estimation  of  ionospheric  delay  errors. 
The  methods  for  determining  ionospheric  delay  errors  for  the  single 
and  dual  frequency  have  been  described  above  in  Section  1.1.2.  2. 
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I.J.3.  5 Other  Parameters 


The  receiver  model  also  consists  of  the  required  interface 
between  itself  and  the  navigation  processor.  It  will  be  necessary  to 
specify  the  output  formats,  resolutions,  time  delays,  frequency  of  mea- 
surements, and  acquisition  times.  These  values  can  be  considered  as 
parameters  of  the  study. 

1.1.4  Reception  Model 

For  the  purpose  of  this  study,  a simplified  reception  model  will 
be  used.  It  will  be  assumed  that  there  is  sufficient  gain  in  the  receiver 
and  user  antenna  to  receive  signals  from  any  satellite  which  is  a fixed 
angle  Y (five  degrees  in  this  study)  above  the  horizontal  phane. 

The  effects  on  satellite  visbility  due  to  platform  pitch 
and  roll  will  be  modelled  as  a decrease  in  the  cone  of  reception  corre- 
sponding to  the  pitch  and  roll  angle.  This  is  a reasonable  model  since 
pitch  and  roll  periods  in  other  simulations  (5)  are  shorter  than  probable 
acquisition  times  (6).  Further  study  will  have  to  be  done  to  see  if  two 
separate  parameters  are  required  to  characterize  the  reduction  of  the 
cone  of  reception.  Two  parameters  may  be  required  since  the  require- 
ments to  re -acquire  a signal  momentarily  lost  may  be  considerably  dif- 
ferent than  the  requirements  to  acquire  a new  signal. 

1.1.5  Integration  of  Other  Sensors 

The  system  model  also  contains  provisions  for  integrating  mea- 
surements from  three  other  sensors  into  the  navigation  algorithm.  These 
sensors  are  the  EM  log,  gyrocompass,  and  Omega  receiver.  The  sys- 
tem and  error  models  for  these  are  taken  from  reference  7.  The  refer- 
ence goes  into  considerable  detail  in  the  development  of  these  models. 
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It  also  cites  the  primary  references  for  the  models  to  be  used  in  this 
study.  These  sensors  are  actually  external  to  the  NAVSTAR/GPS  sys- 
tem. Their  inclusion  is  to  evaluate  the  accuracy  improvement  obtain- 
able from  the  relatively  minimal  cost  of  the  integration  of  these  data 
sources.  It  may  also  be  possible  to  use  the  GPS  data  to  calibrate  the 
other  sensors  for  times  when  GPS  is  not  available. 

1.  1.  5.  1 EM  Log 

The  EM  log  is  an  instrument  which  measures  a ship's  speed 
along  its  longitudinal  axis.  The  measured  speed  is  with  respect  to  the 
water  so  that  ocean  currents  become  a source  of  error. 

1.1. 5. 2 Gyrocompass 

The  gyrocompass  is  an  instrument  aboard  ship  which  indicates 
the  ships  heading.  This  heading  information  then  combined  with  the  EM 
log  data  gives  the  ship's  velocity.  The  pitch  and  roll  information  can 
be  used  to  modify  the  set  of  possible  satellites  which  can  be  tracked. 

1.  1.  5.  3 Omega  Receiver 

Omega  is  a land-based  worldwide  coverage  hyperbolic  radio  navi- 
gation system.  The  Omega  receiver  receives  signals  from  the  transmit- 
ting stations  and  determines  the  user  position.  The  basic  measurable  is 
a phase  difference  from  two  transmitting  stations.  The  system  model 
consists  of  four  Omega  transmitters  so  that  three  independent  lines 
of  position  are  available. 

1.1.6  Scenarios 


The  following  scenarios  have  been  selected  for  use  in  the  analysis 
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of  the  navigation  algorithms.  These  scenarios  were  chosen  at  random 
and  not  contrived  to  emphasize  any  particular  point.  The  scenarios  are 
intended  to  be  two  typical  cases.  They  are  depicted  in  Figures  1-1  and  1-2. 


Test  Scenarios 

#1  Ship  initially  at  45°N  160°W  travelling  due  south  at  10  kts 

at  t = 30  secs,  accelerate  to  30  kts  at  rate  of  . 5 kts /sec 
at  t = 72  secs  turn  starboard  to  360°  at  l°/sec  and  decelerate  to 
15  kts  at  .25  kts /sec 

at  t = 260  secs  turn  port  to  300°  at  3°/sec  followed  by  an  additional 
turn  to  port  of  30°  at  . 5° /sec 
at  t = 350  secs  decelerate  to  stop  at  . 1 kts /sec 

at  t = 530  secs  accelerate  at  . 05  kts /sec  and  turn  port  at  accelerat- 

o 2 

ing  rate  of  . 01  / sec  for  100  secs.  Maintain  turn  rate  at 
l°/sec.  Change  acceleration  to  .2  kts/sec  for  40  secs. 
Maintain  velocity  and  heading  for  40  seconds. 

Total  run:  710  secs 

#2  Aircraft  initially  at  32°N  120°W  travelling  West  at  600  kts. 
Altitude  38,000  ft 

at  t = 0 secs  accelerate  at  2 kts/sec  for  5 min 
at  t = 300  secs  turn  to  port  at  5°/sec  for  36  secs, 
at  t = 360  secs  turn  port  at  . 5°/sec  for  120  secs 
End  at  t = 600  secs 


■) 
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1.2  Error  Model 


# 

i 

The  error  model  ie  a model  of  those  effects  which  tend  to  degrade 
the  accuracy  of  the  NAVSTAR/GPS  navigation  mechanisation.  To  make  a 
meaningful  evaluation  of  any  proposed  navigation  algorithm,  it  is  impor- 
tant to  have  a realistic  error  model  containing  all  of  the  significant  error 
sources.  However,  in  a system  as  complex  as  the  NAVSTAR/GPS  navi- 
gation system,  this  groundrule  would  demand  an  error  model  consisting 
of  hundreds  of  states.  The  error  model  outlined  below  will  attempt  to 
simplify  most  of  the  individual  error  models  to  keep  the  total  number  of 
states  down  while  maintaining  the  integrity  of  the  error  model.  Each 
reduced  state  error  model  will  be  justified  within  the  report  or  the 
references  cited. 

Much  of  the  analytic  work  which  uses  this  error  model  takes  the 
^ form  of  linearized  covariance  analysis.  For  this  reason  it  is  desirable  to 

not  only  make  the  error  models  simple,  but  to  pose  as  many  of  them  as 
possible  as  linear  models;  i.e.,  error  models  whose  behavior  is  described 
by  linear  differential  or  difference  equations.  Models  of  this  form  are  then 
easily  adapted  to  methods  of  linear  analysis. 

1.2.1  Ephemeris  Error  Model 

A complete  error  model  for  satellite  position  and  velocity  uncer- 
tainties is  a large  order  system  containing  many  high- order  gravity  har- 
monic errors.  This  type  of  model  is  used  and  is  necessary  for  accurate 
satellite  position  determination  and  prediction  over  long  time  intervals. 
When  residual  errors  to  accurate  position  predictions  are  considered, 
the  effects  of  the  high-order  gravity  harmonics  are  negligible.  This 
justifies  the  use  of  the  two-body  orbit  equations,  utilising  just  six  states 
per  satellite,  for  propagating  satellite  position  errors.  These  six  states 
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can  be  propagated  from  initial  conditions  by  a state  transition  matrix 
calculated  in  closed  form  (8). 

Propagating  six  states  per  satellite  for  up  to  eleven  visible  satel- 
lites is  still  a large  number  of  states.  Further  reduction  in  the  number 
of  states  is  thus  desirable.  One  way  of  doing  this  is  to  assume  that  the 
position  uncertainty  is  constant  rather  than  growing  with  time.  Since 
the  growth  is  very  small,  choosing  the  maximum  user  equivalent  range 
error  as  the  constant  uncertainty  along  each  axis  will  give  a realistic 
though  somewhat  pessimistic  model  for  the  ephemeris.  Furthermore, 
assume  that  the  errors  are  given  initially  in  the  principal  axes,  so  that 
the  initial  error  covariance  matrix  is  diagonal.  This  can  be  expressed 
by 


0 

O 

O 

Ef  AX]  = 

0 

0 

;E[AX  AXT]  = 

0 (T2  0 

0 0 oZ 

where 

AX  is  the  satellite  position  error  vector 

0 is  the  user  equivalent  range  error  for  Phase  III  specified 
in  Reference  9,  p.  8. 

Another  source  of  error  in  the  ephemeris  model  arises  numeri- 
cally. This  error  is  the  difference  between  the  control  segment  predicted 
satellite  position  and  the  position  computed  by  the  user  from  the  received 
data.  In  reference  (15).  it  is  shown  that  this  error  source  can  be  made 
negligible  by  scaling  of  data  and  careful  algorithm  design. 


A 
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1.2.2  Time  Error  Model 


As  in  the  system  time  model,  the  time  error  model  consists  of 
a satellite  clock  error  model  and  an  atmospheric  delay  error  model.  In 
addition  a user  clock  error  model  must  be  considered.  The  time  error 
model,  like  the  ephemeris  error  model,  is  concerned  with  the  modelling 
of  the  residual  errors  between  actual  time  and  the  system  model  of  time 
(Section  1.  1.2).  Now  the  errors  in  time  and  ephemeris  are  correlated 
inherently  by  the  methods  which  the  control  system  segment  uses  to  es- 
timate the  system  model  parameters  (3).  Ignoring  such  correlations 
may  give  somewhat  optimistic  results,  however  determination  of  such 
correlations  would  require  extensive  simulation  of  the  control  system 
segment  estimation  techniques.  This  is  clearly  not  within  the  scope  of 
the  present  study. 

1.2.2. 1 Satellite  Clock  Error  Model 


The  residual  time  error  consists  of  two  parts,  the  error  in  the 

predicted  correction  parameters  and  the  random  part.  The  errors  in 

the  correction  parameters,  which  represent  bias  error,  frequency  error, 

and  frequency  rate  error,  are  used  to  determine  initial  conditions  for 

the  dynamic  error  model  of  the  satellite  clock.  The  type  of  error  model 

required  for  the  random  part  can  be  determined  from  the  curves  of  the 
* 

Allan  variance  (10)  for  the  clock.  For  the  Phase  III  system.  Cesium 
beam  clocks  will  be  used  in  the  satellites.  These  clocks  will  be  updated 
every  day  so  that  the  error  model  is  based  on  the  twenty-four  hour  vari- 
ance. The  Allan  variance;  of  the  typical  Cesium  beam  clock  data  indicates 

the  fractional  frequency  error  can  be  modelled  as  a white  noise  (i.e.  , 

_ 

The  Allan  variance  is  the  variance  of  the  fractional  frequency  error  as 
a function  of  the  sampling  interval  r.  See  reference  10  for  details. 
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Flicker  noise  and  integrated  white  noise  are  not  important  for  times 
less  than  one  day).  The  satellite  clock  error  model  is  thus 


ft>C50-C) 


u(t) 


where 


and 


is  the  time  error  of  the  nth  satellite  clock 

X2n  is  the  fretluency  bias  error  of  the  n*h  satellite  clock 
0 is  PSD  magnitude  (typically  for  Cesium  Beam  0 = 10 
sec/sec  ) 

u(t)  is  white  noise  with  unity  PSD 


x (0)  is  the  error  in  a 
In  o 

x (0)  is  the  error  in  a, 
cn  1 
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1.2.  2.  2 User  Clock  Error  Model 

The  user  clock  error  model  is  based  on  data  for  typical  crystal 
oscillators.  The  Allan  variance  for  crystal  oscillator  clocks  indicates 
that  the  error  model  over  the  times  of  interest  consists  of  a fractional 
frequency  error  model  with  white  noise,  flicker  noise,  and  integrated 
white  noise.  Following  the  procedures  in  reference  10,  a clock  model 
for  the  user  which  also  includes  a frequency  offset  and  an  aging  coeffi- 
cient is  given  by  the  following  seven- state  description. 
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Tj»  T^,  are  parameters  of  the  Allan  variance  curve  (10) 


Typical  values  for  a crystal  oscillator  are  = . 5 sec, 

4 

T ^ = 80  sec,  T ^ = 5x10  sec 

The  state  is  the  user  clock  error  in  seconds.  The  parameters  T^, 

j t ^ , and  M will  depend  on  the  particular  oscillator  type  chosen  for 
the  user  frequency  standard. 

1.2. 2.  3 Atmospheric  Delay  Error  Model 

The  atmospheric  delay  error  model  consists  of  a residual 
error  model  for  both  ionospheric  and  tropospheric  delays.  The  larg- 
est uncertainties  are  the  ionospheric  model  when  a single  frequency  is 
used.  The  ionospheric  error  model  discussed  here  is  Ifte  residual 
model  for  the  single  frequency  case.  The  dual  frequency  correction 

will  be  considered  to  have  negligible  residual  time  error. 

The  tropospheric  error  model  is  an  uncertainty  in  the  correc- 
tion constant  K (see  Section  1.  1.2.2  above).  The  time  error  in  the 
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incoming  signal  is  a zero  mean  random  constant  AK  times  the  cosecant 
of  the  elevation  angle.  The  variance  of  £K  is  approximately  . IK  (3). 

This  error  source  is  independent  of  the  other  time  models. 

The  ionospheric  residual  error  model  (see  references  11  and  14) 
is  modelled  as  a Gauss-Markov  error  source  which  is  correlated  in  both 
time  and  distance.  The  covariance  of  two  measurements  at  two  different 
times  and  two  different  places  will  be  modelled  as  (12) 

°ZR[r2  * AR,  *R2,t<*>y  *» 

where 


AR2  are  the  prediction  residuals 

AR.  = « esc [Ve^  + (18°)2] 

C is  the  RMS  correction  error 

is  the  elevation  angle  for  the  i**1  measurement  in  degrees 

At  is  the  time  difference  between  measurements 

Ap  is  distance  between  the  ionospheric  "pierce  points" 
of  the  measurements  (The  pierce  point  is  the  point  at 
which  the  transmitted  signal  intersects  the  Ionosphere. 
The  ionosphere  here  is  assumed  to  be  a thin  shell  350  km 
above  the  surface  of  the  Earth) 

F or  the  simulations  in  this  study,  the  functions  1)^  and  can  be 
approximated  over  the  range  of  interest  from  the  data  in  reference  12. 

(The  model  given  in  references  12  and  14  have  two  time  constants,  but 
for  times  in  the  selected  scenarios,  one  is  sufficient. ) 


27 


TJt(T)  = 


e-t/6.  7 hr. 


iw  - .-o'2500  km 


In  order  to  model  this  covariance  in  a linear  system  model,  the 
distance  correlation  will~be  changed  to  a time  function.  This  can  be  done 
since  the  "pierce  point"  of  the  350  km  altitude  ionosphere  shell  moves 
with  a constant  velocity  when  the  satellite  orbits  are  circular.  Thus  re- 
placing p with  /3  • t,  when  p is  constant,  transforms  the  spatial  correla- 
tion into  a time  function.  Combining  TJt  and  ij  then  yields  a Gauss-Markov 
model  with  a modified  time  constant.  So 


R1R2 


AR^R2e 
ARjAR2  e 


- r /6.  7 hr.  - r*  p/2500  km 

-T( — + 2 ) 

T o.  7 hr  2500km’ 

-r/1.  17  hr 


where  <5  Z 2ff  ' (6378  km  + 350  km)/24  hr  S 1761  km/hr.  24  hr  is  the 
apparent  period  of  the  satellite  with  respect  to  an  Earth  fixed  coordinate 
frame.  In  addition,  p can  be  modified  to  account  for  user  motion. 

The  ionospheric  delay  error  for  each  satellite  can  be  modelled 
as  above.  In  addition,  a cross-correlation  of  the  ionospheric  delay 
error  between  the  various  satellites  can  be  computed  using  the  original 
formula  above. 


1.2.3  Receiver  Measurement  Error  Model 

The  exact  form  of  the  receiver  is  as  yet  unresolved.  The  error 
model  of  the  measured  data  is  for  one  channel  of  both  time  of  arrival 
measurements  (for  pseudo-range  or  range  difference)  and  Doppler  mea- 
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surements  (for  pseudo  range  rate).  The  Doppler  errors  are  not 
independent  for  the  case  of  multi-channel  reception.  The  multi-channel 
error  nodel  can  be  put  together  from  copies  of  the  single -channel  model 
with  the  correlations  noted  and  accounted  for. 


1.2.  3.1  TOA  (Time  of  Arrival)  Measurement  Error  Model 


The  primary  error  in  the  TOA  measurement  is  the  user  clock 
error  which  has  been  described  in  Sec.  1.2.  2.  2.  In  addition  to  this  is 
a uniformly  distributed  random  error  corresponding  to  the  measure- 
ment resolution.  This  may  result  from  the  clock  resolution  or  the  finite 
word  length  of  the  time  data.  The  mean  error  will  be  different  for  the 

case  where  the  data  is  rounded  rather  than  truncated.  Let  fit  denote 

res 

the  resolution  error,  then 


E(flt 


res 


(0  for  rounded  error 

. 5 At  for  truncated  error 
res 


E(6t 


res 


res 

12 


where  /St  is  the  smallest  unit  of  time  resolvable.  If  the  resolution 
is  very  fine,  other  measurement  errors  will  limit  how  small  the  TOA 
errors  may  become.  In  either  case  a random  measurement  error  will 
be  included. 


1.2.  3.  2 Doppler  Measurement  Error  Model 

Ideally  a Doppler  measurement  is  an  instantaneous  determination 
of  frequency.  Practically,  this  is  not  possible  so  a count  of  frequency 
over  a short  time  Is  done.  This  can  be  represented  as 
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N 


A 

= $ (“W-  “s,d' 


where 


N is  the  Doppler  count 

CO  .is  the  reference  frequency  CO  = co_  + to 
Rci  Rei  jv  o 

cor  is  the  receiver's  estimate  of  the  transmitted  signal  frequency 


is  an  offset  frequency  so  that  N is  always  > 0 


Wg  is  the  received  frequency  co^  = u)^  + co^ 


is  the  transmitted  frequency 
u>^  is  the  Doppler  shift  frequency 

Now  if  the  time  period  tj  - t^  = 6 t is  short  compared  to  the  domi 
nant  dynamics  of  the  system;  i.e.,  the  user  and  the  satellite,  the  Dopple 
frequency  can  be  assumed  constant  over  fit  and 

CO 

T . 

CO  = — • p 

D C p 

where  p is  the  assumed  constant  range  rate  between  the  user  and  the 
transmitting  satellite,  p = V • u . V is  the  relative  velocity  and  u is 
the  unit  line-of-sight  (LOS)  vector,  and  C is  the  speed  of  light. 

The  Doppler  count  equation  can  then  be  integrated  to  get  N 
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N = (wR  + - coT)8t  - 6t 


The  error  equation  for  AN,  the  Doppler  count  rebidual  is 

x i W_fl  t 

AN  = A6tu>T(l  + £)  - 6tAwT(l  + £)  - A0 


since  p « C,  this  reduces  to 


AN  = Afittoj,^  - flttA^rp  - 6tAA 


The  above  equation  has  assumed  that  time  6t  and  the  frequency 
(WR  + U>q)  are  derived  from  the  same  oscillator  so  that  the  product 
(<i)n  + U)  )fit  is  error  free. 

The  term  A0  is  a function  of  uncertainties  in  user  position,  user 
velocity,  satellite  position,  and  satellite  velocity.  Of  these,  the  satel- 
lite contributions  are  negligible  since  the  ephemeris  data  gives  extreme- 
ly accurate  delta  range  information.  Expanding  tkf)  in  terms  of  the  user 
states  yields  (see  Ref.  7 for  details) 
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where 


u has  been  assumed  to  be  [(X  - X ),  (Y  - Y ),  (Z  - Z )]  • 

u s u s u s 

p is  the  user  to  satellite  range 

X , Y , Z are  the  user  coordinates 
u u u 

X , Y , Z are  the  transmitting  satellite  coordinates 

8 8 S 


m 


Not  included  in  the  above  model  are  the  refractive  errors.  The  refrac- 
tive errors  occur  because  the  transmitting  satellite  is  outside  the  atmo- 
sphere. It  is  shown  in  reference  13  that  the  Doppler  shift  due  to  the  trans- 
mitter and  the  receiver  velocities  are  dependent  on  the  velocity  of  light 
in  their  respective  media.  Ignoring  this  leads  to  a negligible  error  since 
the  light  velocity  difference  is  of  the  order  of  . 03%.  The  other  refrac- 
tion error  is  a change  in  the  direction  of  the  incoming  signal  from  the 
calculated  line  of  sight.  Arguments  similar  to  those  in  Ref.  13  show 
that  this  error  is  negligible  in  the  GPS  geometry. 

Additionally,  there  is  a truncation  or  round-off  error  similar  to 
the  TOA  error  above. 


1.2.4  Other  Sensors 

The  error  models  for  the  EM  log,  gyrocompass  and  Omega  are 
taken  from  Ref.  7. 

1.2.4.  1 EM  Log  Error  Model 

The  state  space  error  model  for  the  EM  log  is 
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where  represents  a random  walk  error  state  and  represents  a 

Gauss-Markov  error  source  with  correlation  time  = . 15  hrs  and  RMS 

2 

value  = . 422  fps  ; u^  and  u^  are  white  noise  with  unity  PSD;  E[xj(0)]  = 
71.23  fps2;  E[x22(0)]  = . 178  fps2. 

The  output  model  is 

y = + x2  + ou3 


where 

a = .472  fps  measurement  noise  (la) 
u^  is  white  noise  with  unity  PSD 

1.2.4. 2 Gyrocompass 

The  state  space  error  model  for  the  gyrocompass  is 

(:;)  ■ C ... 

where 

Xj  is  a random  walk  bias  state  . 013°/^ec(l a) 

x2  is  a Gauss-Markov  error  with  correlation  time  = 23  min 
and  RMS  value  = . 35° 

u j and  u2  ere  white  noises  with  unity  PSD 
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The  output  is  modelled  as 


y = x2  + x2  + <ru3 


where 

!.  15  deg  la  no  maneuvers 
. 52  deg  la  during  maneuvers 
u3  is  white  noise  with  unity  PSD 


1.2.4.  3 Omega  Error  Model 


The  details  of  the  derivation  of  the  Omega  error  model  can  be 
found  in  Ref,  7.  The  state  space  equations  for  the  phase  error  for  each 
station  are: 


where 


Xj  is  a bias  state 


is  a Gaus s -Markov  error 


*3*  f°rm  a ainuaoidal  (periodic)  error  with  period  W 


P = (l/3600)sec 

fa)  = (2ir/12)  hr-1 

2 2 2 * 

O - .00277  centicyclea  /iec 

, 2 . 2 
E[Xj  (0)]  = 8.  5 centicyclea 

r 2 , 2 

Etx^  (0)]  = 5.  centicyclea 

u't)  ia  a white  noiae  with  unity  PSD 

, 2 2 
E[x^  (0)]  = 5.  centicyclea 

E[x^  (0)]  = 1.05  x 10  centicyclea^ 

Converting  thia  phaae  error  into  equivalent  poaition  error  of  the 
uaer  ia  a function  of  the  uaer  location,  Detaila  of  the  calculation  may 
be  found  in  reference  7. 

1.2.5  Platform  Error  Model 

To  account  for  random  motiona  and  accelationa  of  the  uaer,  the 
poaition,  velocity  and  acceleration  atatea  will  be  modelled  aa  random 
w*lke.  The  magnitudea  of  the  random  walk  driving  noiae  variancea  will 
vary  according  to  the  type  of  platform  and  the  acenario. 


*A  centicyclr  ia  1/100  of  a cycle. 
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Section  2.  - SINGLE  FIX  ALGORITHM 


2. 0 Introduction 


This  section  is  concerned  with  the  problem  of  determining  an  ac- 
curate estimate  of  a user  location,  i.  e.  , a navigational  fix,  using  the 
minimum  amount  of  data  required  for  the  fix.  This  first  step  in  the  de- 
velopment of  the  NAVSTAR/GPS  navigation  algorithm  is  called  the  single 
fix  algorithm.  In  some  system  configurations,  the  initial  fix  may  be  made 
8 imply  by  entry  of  present  position  derived  from  other  navigation  aid 
sources  such  as  inertial  navigation  systems,  star  fixes.  Omega,  etc., 
or  from  known  initial  conditions.  The  single  fix  algorithm  is,  however, 
needed  to  establish  initial  conditions  for  some  system  configurations. 
Additionally  it  may  be  used  to  make  the  NAVSTAR/GPS  system  autono- 
mous even  if  other  data  sources  are  available  and  it  may  be  used  in  some 
pseudo-measurement  mechanisations  (e.  g.  the  <*-/3  filter  in  section  3). 

The  single -fix  algorithm  as  presented  in  this  report  consists  of 
two  distinct  parts.  The  first  part  is  the  selection  of  the  available  satel- 
lites to  use  for  the  navigation  fix.  The  second  part  is  the  actual  position 
determination  using  the  time  of  arrival  data  from  the  selected  satellites. 
The  single -fix  algorithms  may  be  mechanized  for  either  position  only  or 
position  plus  user  clock  bias.  The  body  of  the  report  is  concerned  with 
the  latter  case.  The  mechanization  for  the  position  only  algorithm  is 
discussed  in  Appendix  2 of  Reference  16. 

Processing  of  the  measurables  for  the  single-fix  algorithm  falls 
into  two  categories,  batch  and  sequential  processing.  This  section  is 
concerned  with  iterative  techniques  for  batch  processing  of  data.  Se- 
quential processing  mechanizations  fall  more  naturally  in  the  domain 
of  recursive  filtering  algorithms.  As  such  these  techniques  will  be 
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studied  in  subsequent  sections  of  this  report  as  part  of  the  start-up  pro- 
cedures for  the  recursive  filtering  algorithms.  Recursive  filters  which 
can  be  initialized  from  relatively  crude  initial  conditions  may  be  thought 
of  as  having  a built-in  single-fix  capability.  In  this  sense  a separate 
single -fix  algorithm  as  described  herein  may  not  be  needed  to  initialize 
all  of  the  candidate  recursive  filters.  This  does  not,  however,  preclude 
the  use  for  other  purposes. 

This  section  is  organized  in  four  subsections  corresponding  to 
four  areas  of  the  single-fix  algorithm  development  and  analysis.  Section 
2.  1 describes  the  candidate  alert  algorithms,  one  from  the  GDE  proposal  (3) 
and  one  proposed  maximum  volume  using  a known  principle  for  approximate 
minimization  of  GDOP.  Section  2.  2 developes  the  algorithms  for  solving 
the  set  of  non-linear  equations  to  determine  position.  Section  2.  3 describes 
the  techniques  used  to  evaluate  the  algorithms  developed  in  Section  2.  2. 

The  results  from  the  studies  of  Sections  2.  1-2.  3 are  given  in  Section  2.  4. 

All  of  the  computer  programs  used  in  the  analyses  are  given  in  Appendix  1 
of  Reference  16. 

2.  1 Alert  Algorithm  for  Single -fix 

The  measure  of  a good  alert  algorithm  is  the  value  of  GDOP  (6) 
for  the  satellites  which  it  selects.  The  nature  of  GDOP  does  not  allow 
for  any  exact  methods  of  GDOP  minimization  short  of  an  exhaustive 
enumeration  of  all  possibilities.  With  from  six  to  eleven  satellites 
visible  at  any  point  in  time,  an  exhaustive  enumeration  may  take  too 
much  computation  time.  For  this  reason,  two  approximate  GDOP  mini- 
misations are  also  being  considered  as  candidates  for  the  alert  calcula- 
tions. The  first  technique  is  the  one  developed  by  GDE  for  its  proposal  (3). 
The  second  technique  is  based  on  an  approximate  maximization  of  the 
tetrahedron  enclosed  by  four  satellites.  The  form  of  the  GDOP  calcula- 
tions used  will  be  described  in  Section  2.  3. 


2.  1.  1 GDE  Alert  Algorithm  (3,  p.  1-39  to  1-41) 


This  algorithm  gives  criteria  for  the  selection  of  the  four  satel- 
lites to  be  used.  The  minimizations /maximizations  to  be  performed 
involve  only  dot  products  of  vectors  and  so  are  simple  and  fast  compu- 
tationally. This  method  proceeds  as  follows: 

Satellite  No.  1 selection:  Choose  the  satellite  which  is  closest 
to  the  zenith.  This  is  found  by  maximizing  the  dot  product  of  the  unit 
vector  to  the  user  with  the  unit  vector  from  the  user  to  the  satellites, 
st 

1 Satellite  No.  = i such  that  U • U * U * U 

“ i ~u  — j -u 

where 

Uu  is  the  ECI  unit  vector  to  the  user 

U,  is  the  unit  vector  along  the  line  of  sight  from  the  user  to 
th 

the  k satellite  ( k=i,  j ) 

j varies  over  the  visible  satellites 

Satellite  No.  2 selection:  The  second  satellite  is  chosen  as  the 
satellite  among  those  visible  which  is  closest  to  the  horizon  but  not 
in  the  same  orbit  plane  as  the  first  selected  satellite. 

2n<*  Satellite  No.  ■ i such  that 

1)  Uj  * with  j varying  over  the  visible  satellites 

2)  i is  not  in  the  same  orbit  plane  as  the  first  selected  satellite. 

Satellite  No.  3 selection:  The  third  satellite  is  selected  as  the 
one  among  the  visible  satellites  which  is  closest  to  being  orthogonal  to 
both  the  first  and  second  satellites  selected. 

3r<*  Satellite  No.  = i such  that 

|v  <^i  *M*|v  (-i  * -2>l 
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> where 

j varies  over  the  visible  satellite  and  is  the  unit  vector 

from  the  user  to  the  satellite  chosen  as  Satellite  No.  1 and 

U,  is  the  unit  vector  from  the  user  to  the  satellite  chosen  as 
* “ 2 

Satellite  No.  2. 

Satellite  No.  4 selection:  This  satellite  is  chosen  as  the  one 
which  is  closest  to  the  vector  sum  of  the  first  three. 

4**1  Satellite  No.  = i such  that 

Hr  <y,  + h2  + y3)  * Hj  • <Hi+y2  + y3) 

j varies  over  the  visible  satellites 

where 

is  the  unit  vector  from  the  user  to  the  satellite  chosen  as 
Satellite  No.  3. 

It  should  be  noted  that  the  reference  contains  some  obvious  errors  which 
have  been  corrected  in  this  summary.  The  above  described  algorithm 
is  the  one  which  was  used  for  purposes  of  comparison. 


2.1.2  Maximum  Volume  Alert  Algorithm 


where 


j varies  over  the  visible  satellites. 

Satellites  No.  2,  3,  and  4 selection:  These  satellites  are  chosen 
as  a group  by  iterating  once  through  all  of  the  remaining  satellites.  The 
procedure  is  as  follows. 


A)  Pick  three  visible  satellites  and  compute  the  tetrahedral 
volume  (see  Figure  "Hexahedron  of  satellite  and  user". 

B)  Pick  another  visible  satellite,  if  there  are  no  more,  then 
procedure  is  finished. 

C)  Compute  the  volume  of  the  tetrahedrons  with  the  new  satel 
lite  successively  replacing  satellites  No.  2,  3,  and  4. 

Pick  the  configuration  with  the  largest  volume  and  label 
those  included  satellites  as  the  new  No.  2,  3,  and  4. 

(Go  to  Step  B). 

This  algorithm  for  finding  the  best  satellites  creates  a monotonically 
decreasing  GDOP.  This  procedure  converges  to  a local  minimum  and 
it  is  certainly  dependent  on  the  order  of  the  satellites  chosen.  In  the 
analysis  section,  the  order  taken  is  arbitrary  so  that  the  analysis  repre- 
sents a lower  bound  on  the  ability  of  the  technique  to  choose  a good  con- 
stellation (i.  e. , an  upper  bound  for  an  achievable  GDOP). 

The  volume  computation  is  a simple  calculation.  The  volume  may 
be  computed  using  the  scalar  triple  product. 


* »/14)  l»,2  ' (»,j  « »,4)l 
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where 


Vj  is  the  volume  of  the  tetrahedron  labeled  I (See  Figure  2.  1) 

a^  is  the  vector  from  Satellite  No.  1 to  Satellite  No.  j. 

] | indicates  absolute  value 

With  this  algorithm  at  most  twenty -two  volume  computations  need  be 
done  when  eleven  satellites  are  visible. 

2.  2 Single  Fix  Algorithm 

This  section  describes  the  candidate  solution  techniques  for  the 
simultaneous  non-linear  equations  which  desctibe  the  user  location. 

The  candidate  algorithms  presented  here  will  be  evaluated  on  the  basis 
of  convergence,  sensitivity  to  initial  conditions,  amount  of  computation 
required,  and  accuracy  in  the  presence  of  noisy  inputs.  The  evaluation 
techniques  and  results  will  be  described  in  the  next  two  sections.  All 
of  the  analysis  in  this  and  subsequent  sections  is  for  position  plus  user  bias 
fixes.  The  formulation  for  the  position  only  fixes  is  contained  in  Appen- 
dix 2 of  Reference  16. 

To  compute  a fix  from  satellite  range  data,  measurements  from 
Uiree  satellites  (for  a two-dimensional  fix)  or  four  satellites  (for  a three- 
dimensional  fix)  are  required.  For  this  study  it  has  been  assumed  that 
all  fixes  are  three  dimensional.  The  basic  measurement  to  be  used  for 
position  computation  is  the  signal  time -of-ar rival  (TOA).  To  use  a TOA 
in  a navigation  algorithm,  the  following  pieces  of  information  are  required 
for  each  measurement: 

A)  Position  of  transmitting  satellite  at  time  of  signal  transmission 

B)  Time  of  transmission  of  the  received  signal 

C)  Estimate  of  the  deterministic  time  delays 
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With  this  information,  the  user  may  use  the  TOA's  from  selected  satel- 
lites  to  determine  its  position.  The  equations  to  be  solved  can  be  devel- 
oped as  follows: 

The  measured  time  of  arrival  from  satellite  i,  AT . , is  given  by 
ATj  = (1/c)  ||  ^-X  ||  +Ti  + t.+b  + w 


where 

is  the  position  of  the  i**1  satellite  at  the  time  of  transmission 

"1 

X is  the  user  position 

T.  is  the  transmission  time 
1 

ti  is  the  deterministic  delay 
b is  the  user  bias  (clock  and  electronic  delay) 
w is  the  measurement  error  due  to  receiver  error,  random 
atmospheric  delays,  ephemeris  errors,  satellite  electronic 
random  delays,  etc. 
c is  the  speed  of  light 

When  a set  of  AT  have  been  measured,  a position  fix  can  be 
obtained.  Two  cases  must  be  considered  for  collecting  sets  of  data. 

The  two  cases  correspond  to  the  single  channel  and  four  channel  re- 
ceiver configurations.  In  the  first  case  there  will  be  motion  of  the  user 
platform  between  measurements;  in  the  second  case,  there  is  no  motion. 
The  user  motion  in  the  first  case  will  lead  to  some  error  unless  an  ex- 
ternal velocity  determination  is  available. 

The  set  of  measurements  leads  to  the  simultaneous  non-linear 

e 

equations.  The  candidate  solution  techniques  are  presented  below.  They 
are  all  Iterative  techniques. 


t 
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2.  2.  1 Newton- Raphs on  Method 


The  Newton -Raphs on  method  is  a well  known  procedure  for 
solving  sets  of  simultaneous  non-linear  equations  (17).  The  system 


of  equations  to  be  solved 

here  is: 

fjQC,  b)  = (1/c) 

I!  Sj  - x 

II 

il 

+ b + (Tj  - ATj) 

= 0 

f2(X,  b)  = (1/c) 

I'  52-2 

1! 

Il 

+ b + (t2  - at2) 

= 0 

f3(X,  b)  = (1/c) 

!!  s3  - x 

If 

ii 

+ b + (t3  - at3) 

= 0 

f4(X-  b)  = (1/c) 

!!  s.  -x 

—4  — 

U 

ii 

+ b + (t4  - at4) 

= 0 

Note  that  if  there  is  a bias  in  the  system  such  that  all  of  the  AT  are 
off  by  a fixed  constant  increment,  this  looks  like  a user  bias  b.  The 
user  and  system  biases  are  inseparable  and  only  {he  sum  may  be  esti- 
mated for  navigational  purposes.  Lumping  the  two  biases  into  b creates 
no  problem.  The  above  set  of  equations  will  be  used  throughout  this 
section  as  the  set  of  equations  to  be  solved. 

The  Newton-Raphson  method  proceeds  from  one  estimate  to 
the  next  according  to 


Xn+1 

= 

X 

n 

+ k 

n 

yn+l 

= 

yn 

+ 1 

n 

z . . 

- 

z 

+ m 

n+l 

• 

n 

n 

bn+l 

s 

b 

n 

+ P 

*ii 

where 

*i»  y are  the  user's  position  estimate  at  the  end  of  the  ;th 
iteration  and  bj  is  the  estimate  of  the  user  bias  at  the  ith  iteration. 
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where 


y denotes  the  gradient  of  the  function  evaluated  at  the  point 
indicated  by  the  argument. 

These  equations  can  be  solved  by  a Gaussian  elimination  scheme. 
No  matrix  inversion  is  required.  Other  than  add,  subtract,  multiply, 
and  divide,  the  only  computational  operation  required  is  a square  root. 

In  the  GPS  geometry,  the  Newton -Raphs on  method  converges  very  well 
with  large  initial  condition  errors. 

2.  2.  2 Non-linear  Gauss-Seidel  Iteration 


The  basic  Gauss-Seidel  iteration  philosophy  is  easily  applied 
to  sets  of  non-linear  equations.  The  simplicity  of  the  iteration  pro- 
cedure makes  this  algorithm  a candidate  for  use  in  the  single -fix  com- 
putation. 

The  Gauss-Seidel  approach  to  sets  of  simultaneous  non-linear 
equation  is  as  follows: 

Let  f.  (xi , x , . . . , x ) - 0 i = i , 2,  . . . , n be  a set  of  si- 

114  n th 

multaneous  non-linear  equations.  During  the  k iteration,  the  update 

of  x.  comes  from  the  i**1  equation  by  solving  for  x^  in 


f.  (X, 


Oc)  _ (k) 


(k) 


Xi+1 


(k-1) 


where  the  superscript  (j)  denotes  the  value  of  the  variable  from  the  j**1 
iteration. 

For  this  problem,  the  iteration  is  obtained  from  the  set  of  equa- 
tions for  the  navigation  problem  (see  Section  3.1).  The  k^  iteration 
is  given  by: 
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x * x,1illAT1.I1.bk.1W2.(yil-Tl./.(.ll-%.tlV/2 

= y»2 -(<AT1  "T1 'bk-l)c,Z  ' (x«2'Xk,Z  ’ (‘.2"*k-l,2,1/2 
% ■ %3i((AT3-T3-bk-l>c)2-«‘.3-V2-'y.3-yk,2»1/2 

bk  ■ AT4-T4-i/c((,t4-xk,2  + ,y.4-yk)2  + (x.4-V2l1/2 

where  x .,  y .»  and  z . are  the  ECI  coordinates  of  the  i^1  satellite  being 
8*81  81 

used. 

This  technique  is  particularly  simple  since  there  are  no  simul- 
taneous equations  to  solve.  There  is,  however,  some  additional  logic 
to  resolve  the  sign  in  the  first  three  equations. 


2.2.3  Successive  Linearizations  of  Measurement  Matrix 

i 

This  method  is  similar  to  the  one  in  (18)  except  that  the  same 
data  will  be  used  in  the  iteration  instead  of  new  data  points.  Because 
of  this,  there  will  actually  be  two  simultaneous  iterations  in  progress. 
Some  efficient  means  for  calculating  initial  conditions  for  the  iteration 
must  be  developed.  The  basic  equations  for  this  method  are  the  linearized 
perturbed  measurement  equations,  i.  e. , 


AT, 


in  . bfig-b» 


x=x 


Ax 


Bfj  (X,b) 
~ 


y*y. 


Ay  + 


(X,  b) 
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z=z 


Az 


6ft(X,  b) 
6b 


b=b 


Ab  , i * 1,  2,  3,  4 


n 
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where  A refers  to  transit  time  of  the  satellite  signal  from  satel- 

lite i.  This  can  be  expressed  as 

[ATt]  = F(Xn.  bj  [AX.  Ab]T 

Then 

[AX.  Ab]T  = F_1(Xn*  V [ATt] 

where  [AX.  Ab]  is  the  update  vector. 

Given  the  initial  conditions  of  the  iteration  viz.  [X  , b ] and 

— o o 

F (X  , b ) the  iteration  proceeds  as  follows. 

— o o 

1.  Compute  vector  [AT,^]  where 

[ATt]  = ATt  - Tj  - [(1/c)  1!  Sj  - xn  !!  n + b] 

2 . Compute  F 

3.  Compute  new  F * 

F * = F (21  - F F”\)  where  I is  the  identity 
n n-1  n-1 

4.  Compute  [AX  , b ] 

~n  n 

[AXn,Abn]T  . f;1  fiTTl 

5.  Compute  [X  , b ] 

n n 

t2„-  b„]  • »„.!)*(“„•  *bn] 

6.  Check  convergence  criteria  and  end  iteration  or  go  to  1. 
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2.  3 Algorithm  Analysis 


Three  areas  of  algorithm  analysis  are  addressed  in  this  section. 

The  first  is  the  analysis  of  the  two  candidate  alert  algorithms.  The  cri- 
teria for  selection  of  one  algorithm  over  the  other  include  lower  average 
GDOP,  computational  complexity,  and  length  of  time  to  do  the  computa- 
tion. The  second  area  of  analysis  is  the  convergence  of  the  numerical 
techniques  presented  in  Sections  2.  2.  1-2.  2.  3.  The  analysis  will  attempt 
to  determine  regions  of  convergence  for  each  algorithm.  The  third  area 
of  analysis  is  the  effect  of  noisy  measurements  on  the  single  fix  algorithms. 
The  possible  effects  are  accuracy  of  the  resulting  fix  and  a possible  change 
in  the  convergence  properties.  The  analysis  techniques  are  described  in 
this  section  with  the  results  in  Section  2.  4. 

2.  3.  1 Alert  Algorithm  Analysis 

The  relative  measure  of  how  accurately  position  may  be  deter- 
mined using  noisy  measurements  from  different  satellite  constellations 
is  called  GDOP  (Geometric  Dillution  of  Precision).  This  measure  is 
a static  quantity  which  is  valid  only  for  a single  set  of  measurements. 

It  represents  an  error  multiplication  factor  relating  the  uncertainty  in 
the  measured  TOA  to  the  resulting  least  squares  position  determination. 

A complete  derivation  of  GDOP  can  be  found  in  Reference  (6).  The  fol- 
lowing formula  for  computing  GDOP  is  derived  in  Reference  (6). 

T -11/2 

GDOP  = (Trace  (F  F)  V 

where 


F is  the  matrix  of  partial  derivatives  («ee  Section  2.  2.  3 above) 
evaluated  at  the  user's  actual  location. 
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The  analysis  of  the  alert  algorithms  uses  the  satellites  in  their  nominal 
orbits. 

To  help  determine  which  alert  algorithm  is  to  be  used,  an  average 
GDOP  value  over  a grid  will  be  computed  using  each  alert  algorithm. 

The  grid  selected  takes  advantage  of  the  symmetries  of  the  problem,  so 
that  only  one -sixth  of  the  Earth  is  considered.  The  grid*  is  every  10° 
in  latitude  from  the  equator  to  the  North  Pole,  every  10°  in  longitude  for 
a total  of  120°,  and  in  time,  every  15  minutes  for  one  and  a half  hours. 
The  average  value  will  be  computed  for  the  GDE  alert  and  the  Max.  Vol- 
ume alert.  The  optimum  achievable  GDOP  for  a given  constellation  can 
be  computed  through  an  exhaustive  search, 

2.  3.  2 Single  Fix  Convergence  Analysis 

Each  of  the  single -fix  algorithm  solution  techniques  has  been  pro- 
grammed for  checkout  purposes.  To  check  the  convergence  of  the  algo- 
rithms, several  sets  of  initial  conditions  were  used  with  varying  error 
magnitudes.  The  convergence  checks  were  also  made  over  varying 
geometric  conditions,  i.  e. , different  user  locations  and  different  times 
to  include  the  changing  satellite  positions. 

2.  3.  3 Effects  of  Noisy  Measurements 

Since  the  measurements  obtained  by  any  GPS  receiver  are  noisy, 
that  is  the  measured  time  is  not  an  exact  indication  of  the  signal  transit 
time,  the  single-fix  algorithms  should  be  checked  for  accuracy  and  con- 
vergence in  the  presence  of  noise.  For  the  analysis  of  the  single-fix 
algorithm,  the  measurement  noise  will  be  modelled  as  a Gaussian 

* 

The  grid  consists  of  654  points. 
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distributed  random  variable  added  to  the  simulated  signal  transit  time. 
The  noise  can  be  assumed  to  be  zero  mean  with  variance  specified  by 
the  phase  III  specification  (6).  The  theoretical  lower  limit  on  the  navi- 
gational error  variance  is  the  measurement  noise  variance  times  the 
GDOP  factor  squared.  This  limit  will  be  reached  if  the  numerical  algo- 
rithms do  not  introduce  errors  which  are  comparable  in  magnitude  to 
this  basic  limitation. 

To  investigate  the  effects  of  noisy  measurements,  a Monte  Carlo 
type  checkout  will  be  done.  The  procedure  will  be  to  select  a user  loca- 
tion, then  simulate  the  measureables.  The  simulated  measurables  are 
given  by 

Tj  = !!  - X 1!  + b + w 

where 

Sj  is  the  simulated  actual  location  of  the  i**1  satellite 

X is  the  simulated  actual  location  of  the  user 

b is  the  simulated  actual  bias 

w is  a random  number  which  is  Gaussian  distributed  with  zero 
mean  and  variance  equal  to  phase  III  specification  for  system 
error  variance,  w is  generated  by  subroutine  GAUSS. 


At  each  user  location,  about  a hundred  fixes  will  be  made  with 
simulated  measurables  as  above.  The  error  will  then  be  computed  as: 


N 


f . il  J v H.l/2 

TIMS  ( N i«i  — -2^") 
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where 


N is  the  number  of  fixes  using  simulated  noisy  data 
X is  the  simulated  actual  user  location 

X.  is  the  computed  fix  using  the  i**1  set  of  noise  measurables 
CRMS  root  square  error 


The  variance  of  the  noise  will  also  be  increased  to  account  for 
inaccuracies  in  the  receiver  measurement.  Good  single  fix  algorithms 
should  not  be  sensitive  to  the  magnitude  of  the  noise  as  long  as  the  mag 
nitude  is  not  unreasonable.  Results  of  this  analysis  are  in  section  2.  4. 


2.  4 Results  of  the  Analyses 

In  this  section,  the  numerical  results  of  the  analyses  described 
in  Sections  2.  3.  1-2.  3.  3.  AH  of  the  computer  programs  to  obtain  the 
numerical  results  are  given  in  Appendix  1 of  Reference  (16). 

2.4.1  Alert  Algorithm  Results 

The  results  of  computing  the  average  GDOP  using  each  of  the 
alert  algorithms  with  a 5°  elevation  horizon  described  in  Section  2.  1 
is  given  in  Table  2.  4. 1.  The  average  was  taken  over  the  space-time 
grid  described  in  Section  2.  3.  1. 


Table  2.  4.  1 Alert  Algorithm  Comparison  Results 


Algorithm 

GDE  Proposal  Alg 

Maximum 

Volume 

Optimum  Configuration 

Average 

GDOP 

47.6* 

2.  84 

2.73 

Max 

Deviation 

3088* 

1.67 

— 

Excluding  singular  points  of  algorithm 
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A more  detailed  look  at  the  GDOP  analysis  is  presented  in  the 
histograms  of  Figures  2.  2 - 2.4.  Using  the  grid  describedin  Section 

4-'  ’ 

2.  3.  1,  Figure  2.  2 presents  a histogram  of  the  GDOP  deviations  of  the 
GDE  proposal  alert  algorithm  constellation  selection  from  the  optimal 
constellation  selected  by  exhaustive  enumeration.  Figure  2.  3 does  the 
same  for  the  maximum  volume  algorithm.  The  distribution  of  the  GDOP 
for  the  optimally  selected  constellations  is  presented  in  Figure  2.  4.  The 
histograms  show  that  the  GDE  proposal  algorithm  selects  poor  constel- 
lations (i.  e.  , GDOP  deviation  > 7)  in  what  is  felt  to  be  an  unacceptably 
large  percentage  of  the  cases.  In  addition,  the  GDE  proposal  algorithm 
produced  some  constellations  for  which  GDOP  did  not  exist,  (i.  e.  , 

GDOP  = •).  This  situation  may  result  since  the  satellites  are  in  nominal 
orbits  with  perfect  symmetry.  Although  in  practice  this  condition  prob- 
ably would  not  arise,  its  possibility  is  a shortcoming  of  the  algorithm. 

To  compare  the  computation  time  of  the  Maximum  Volume  algo- 
rithm with  the  optimum  selection  based  on  volume  maximization,  the 
number  of  volume  computations  should  be  noted.  For  the  Maximum 
Volume  algorithm  the  number  is  3*(N  - 4)  + 1 where  N is  the  total 
number  of  visible  satellites.  For  the  optimum  exhaustive  enumeration, 
the  number  can  be  computed  from  the  binomial  coefficient.  Table  2.  4.  2 
compares  the  numbers  for  the  minimum  number  of  satellites  visible,  6, 
the  typical  values,  8 and  9,  and  the  maximum  number  11.  The  amount 
of  logic  and  hence  the  amount  of  computation  per  satellite  set  required  for 
the  two  mechanizations  is  not  significantly  different. 
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Figure  2.  2 Histogram  of  GDOP  Deviations  for  GDE  Proposal  Algorithm. 
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Table  2.  4.  2 Volume  Computations  Comparison 


No.  of  Visible 
Satellites 

Volume  Computations 
Maximum  Volume 
Algorithm 

Volume  Computations 
Optimum 

6 

7 

15 

8 

13 

70 

9 

16 

126 

11 

22 

330 

2.4.2  Convergence  Analysis  Results 

The  convergence  of  the  three  single  fix  algorithms  was  checked 
for  several  user  locations  and  initial  conditions.  Early  in  the  analysis 
it  was  apparent  that  the  non-linear  Gauss -Seidel  iteration  was  not  well 
suited  for  the  single  fix  algorithm.  The  convergence  was  found  to  be 
very  sensitive  to  initial  conditions,  so  the  method  was  not  considered 
for  further  analyses. 

The  computer  program  written  to  run  convergence  test  cases 
is  called  ITERTEST.  The  program  is  documented  in  Appendix  1 of 
Reference  (16).  This  is  the  same  program  which  will  be  used  to  study 
the  effects  of  noisy  measurements.  The  convergence  checks  are  made 
by  setting  the  noise  variance  to  zero  and  the  number  of  iterations  to  one. 

Table  2.4.  3,  Convergence  Analysis  Sample  Results,  contains 
some  sample  results.  The  initial  conditions  for  the  iterations  were 
chosen  to  correspond  to  an  octant  of  uncertainty.  The  center  of  the 
Earth  was  chosen  in  an  attempt  to  arrive  at  a universally  acceptable 
starting  point. 
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Table  2.  4.  3 Convergence  Analysis  Sample  Results. 


Run 

User  Location 
(Lat.  Long. ) 

Iteration  Initial 

Converged? 

Condition  (in 
Lat.  Long,  or 
Center  of  Earth 
Indication 

NR 

Iteration 

Meas.  Matrix 
Iteration 

Yes 

No 

Yes 

No 

1 

45°N 

160°W 

Center  of  Earth 

X 

X 

2 

45°N 

160°W 

0°N  1 80°W 

X 

X 

3 

45°N 

160°W 

90°N 

X 

X 

4 

45°N 

160°W 

0°N  90°W 

X 

X 

5 

32°N 

120°W 

Center  of  Earth 

X 

X 

6 

32°N 

120°W 

90°N 

X 

X 

7 

32°N 

120°W 

0°N  120°W 

X 

X 

8 

32°N 

120°W 

0°N  90°W 

X 

X 

9 

89°S 

20°W 

Center  of  Earth 

X 

X 

10 

89°S 

20°W 

90°S 

X 

X 

11 

89°S 

20°W 

0°N  0°W 

X 

X 

12 

89°S 

20°W 

0°N  90°W 

X 

X 

2.4.3  Noisy  Measurement  Analysis  Results 

The  Newton -Raphson  iteration  and  the  successive  linearization 
iteration  methods  were  analyzed  with  noisy  inputs.  The  analysis  was 
done  with  main  program  1TERTEST.  Results  are  presented  in  Table 
2.4.  4,  Measurement  Noise  Analysis  Sample  Runs.  Initial  conditions 
were  chosen  to  correspond  to  cases  where  the  iterations  both  converged. 

If  there  is  no  error  or  negligible  error  introduced,  the  RMS 
errors  will  be  GDOP  times  the  input  noise  sigma.  For  the  purpose  of 
this  analysis,  only  one  hundred  iterations  were  done  to  save  computer 
costs.  The  true  error  is  within  twenty-five  percent  of  the  indicated 
error  at  a ninety  percent  confidence  level  when  one  hundred  samples 
are  used.  (X^  test) 


58 


> Conclusions  and  Recommendations. 

The  following  conclusions  and  recommendations  have  been  drawn 
as  a result  of  this  portion  of  the  study. 

A.  The  GDE  proposal  alert  algorithm  is  not  acceptable. 

It  is  recommended  that  an  algorithm  based  on  volume 
maximixation  be  used.  Depending  on  how  much  time  can 
be  allocated  to  this  task  either  the  approximate  or  the 
exhaustive  algorithm  should  be  used. 

B.  The  convergence  of  the  iterative  algorithm  is  not  af- 
fected by  noisy  inputs.  Furthermore,  they  appear  to 
converge  with  a close  approximation  to  the  accuracy  of 
a least  squares  solution. 

C.  If  it  is  desirable  to  have  a single -fix  algorithm  which 
converges  from  an  initial  condition  of  an  octant  of  the 
Earth  uncertainty,  the  successive  linearizations  of  the 
measurement  matrix  method  (sec.  2.  2.  3)  should  be  used 
with  initial  conditions  of  the  center  of  the  Earth. 
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Section  3.  - FILTER  DEVELOPMENT 
3. 0 Introduction 

The  purpose  of  this  section  is  to  present  the  results  of  the  filter 
algorithm  development  task  of  the  study.  There  were  two  objectives  in 
this  task.  The  first  was  to  establish  the  optimum  filter  for  the  NAVSTAR/ 
GPS  receiver  based  on  the  models  in  Section  1.  The  second  objective 
was  to  outline  the  form  of  several  candidate  sub-optimum  filters.  The 
details  of  the  sub-optimum  filters  are  not  established  until  the  end  of 
the  covariance  analysis  section. 

The  optimum  filter  will  include  all  of  the  significant  error  sources 
identified  in  the  modeling  section  of  this  report.  This  will  lead  to  a filter 
with  a very  large  number  of  states,  in  fact  too  large  to  be  considered  for 
actual  implementation  in  an  operational  receiver.  The  purpose  of  estab- 
lishing this  optimum  filter  is  to  set  a bound  on  the  obtainable  accuracy. 

In  addition,  this  optimum  filter  will  become  the  reference  system  for  use 
in  the  covariance  analysis  of  the  next  section. 

Several  sub-optimum  filters  will  be  presented  which  vary  greatly 
in  both  complexity  and  computational  burden.  Each  of  the  filters  is  re- 
cursive in  nature  to  take  advantage  of  the  large  amount  of  data  available. 

For  the  purpose  of  this  section,  only  the  form  of  the  various  filters  can 
be  considered.  Determination  of  the  state  vectors  can  only  be  done  through 
sensitivity  analysis.  Thus  an  evaluation  of  both  accuracy  and  computational 
burden  will  be  presented  in. -a  subsequent  section. 

The  problem  of  actual  operational  mechanisation  for  a given  filter 
will  not  be  discussed  here.  By  this  it  is  meant  that  for  the  purpose  of 
this  report,  for  example,  no  distinction  will  be  made  between  a standard 
or  square  root  formulation  for  the  Kalman  filter.  No  mention  will  be  made 
of  sequential  vs.  batch  processing  except  in  the  cases  where  the  nature  of 
the  filter  demands  one  type  or  the  other. 
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3.  1 Optimum  Filter /R*»f«n»«tf><»e  System 

This  section  will  describe  in  detail  the  optimum  filter.  The  de- 
scription will  include  a summary  of  the  Kalman  filter  equations  and  a 
complete  description  of  the  models  used.  The  model  description  is  a compen- 
dium of  the  models  in  Section  1.  It  serves  a definite  purpose  of  its  own 
since  this  is  the  first  time  that  all  of  the  models  are  brought  together  and 
the  interrelationships  shown.  For  an  understanding  of  the  individual  mod- 
els, it  is  required  to  consult  Section  1. 

The  filters  to  be  described  in  this  and  following  sections  are  bas- 
ed on  linearized  equations.  The  state  equations  for  the  filters  can  be  ex- 
pressed by  linear  relationships.  The  measurements  are,  however,  highly 
nonlinear  functions  of  both  time  and  user  location.  Thus  in  order  to  apply 
the  results  from  linear  system  theory  to  the  estimation  problem  at  hand, 
the  measurement  equations  must  be  linearized  about  some  point.  In  an 
actual  operational  system,  this  point  can  only  be  the  current  best  estimate 
of  the  user  location.  The  procedure  for  the  filter  implementation  is  as 
follows: 

1)  Using  the  current  best  estimate  of  the  state  (either  from  init- 
ial conditions  or  the  value  extrapolated  from  the  previous  estimate) 
and  the  nonlinear  measurement  equations,  compute  the  expected 
values  of  the  measureables. 

-k  = - (<Pk,k-l2k-l* 

2)  Difference  this  expected  measureable  from  the  actual  measured 
value.  This  resulting  value  is  the  measureable  for  the  linearized 
(or  error  state)  equation 

^-*k  * *k  - K 

O 
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3)  Use  the  difference  computed  above  as  the  measurement  in  the 
linearized  filter  equations  to  estimate  the  error  states,  i.  e.  the 
difference  between  the  true  value  and  the  current  estimate. 


4)  Add  the  estimated  error  states  to  the  current  estimate,  X*  , to 
produce  a new  best  estimate.  This  is  sometimes  referred  to  as 
resetting  the  state. 


= ^k.k-l  ^k-1  + -k 

The  linear  systems  descriptions  which  follow  are  the  system 
equations  for  the  errors  of  the  best  estimate.  This  system  is  called 
the  error  system. 

3.1.1  Optimum  Filter  Equations  f 1 91 

The  optimum  filter  equations  are  optimum  in  the  sense  that  they 
provide  the  minimum  variance  estimate  for  a linear  system.  The  opti- 
mum (or  Kalman)  filter  equations  can  be  written  for  either  continuous  or 
d .icrete  systems.  The  nature  of  the  data  from  the  NAVSTAR/GPS  re- 
ceiver dictates  a discrete -time  measurement  and  estimate  update.  On 
the  other  hand,  the  system  has  a continuous  time  system  description  from 
which  a transition  matrix  may  be  computed.  The  system  model  will  be 
presented  as  a continuous  system  for  purposes  of  exposition.  The  filter 
equations  will  be  presented  in  this  and  subsequent  sections  as  discrete 
time  systems  where  the  transition  matrices  r re  computed  using  the  con- 
tinuous time  representation  given. 

The  linearised  differential  equation  for  the  reference  system  is 


identical  to  the  error  system  given  by 


x = Fx  + Gu 


where 


(3.  1) 


x is  the  error  system  state  vector 
F is  the  reference  system  matrix 
G is  the  reference  system  input  distribution 
u is  a vector  of  independent  gaussian  white  noise 
inputs  each  with  unity  PSD. 

The  vector  x has  dimension  n and  u has  dimension  m.  The 
matrices  F and  G are  dimensioned  conformably.  The  linearized  ob- 
servation process  is 

A-k  = Hk*k  + % (3.2) 

where 

A 

^ i®  the  vector  of  predicted  observations 
*®  ike  vector  of  actual  observations 
i#  the  vector  of  observations  for  the  error  system 
H^is  the  error  system  measurement  matrix 
V a vector  of  independent  gaussian  white  noise 
measurement  errors. 

In  addition  z and  v have  dimension  r with  H dimensioned 
conformably.  The  usual  assumptions  on  the  random  processes  are  made 

and  given  by  _ r „ , r 

y E[u(t)]  = E[v(t)]  = 0 

Efe(t)uT(T)]  = 

F[v(t)vT(T)]  = 

E[u(t)vT(T)]  = 


I6(t  - T) 

R(t)fi(t-T)  (3.  3) 

0 
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E[x(t0)uT(t)  ] = Efi^o)  vT(t)  ] = 0 t * tQ 
E[x(t0)xT(t0)]  = P(tQ)  * 0. 

The  state  transition  matrix  <p  for  the  constant  coefficient  refer- 
ence system  of  Eq.  (1)  is  defined  by 


$(t)  = F<p(t) 

(3.4) 

or  solving  for  the  time  invariant  case 

<P(t)  = exp  (Ft) 

(3.  5) 

With  these  definitions,  the  Kalman  filter  equations  for  the 
variance  estimate  x of  the  error  system  state  vector  x , 

minimum 
at  the  k1*1 

measurement  time  are 

K ■ ■ vi*  - 

(3.6a) 

and 

*k  “ ®k, k>lPk- A. k-l+  °k 

(3.  6b) 

“k * pX(HkpkHk  * V‘‘ 

(3.6c) 

pk  ■ t‘-KkHkiI'ki-KkHk>T  + KkRkKk 

(3.  6d) 

At 

Qk  = \ v>(At  - r ) G G<pt(  At  - r)  d r 

(3.  6e) 

0 


where  At  is  the  time  between  measurements  k and  k - 1. 

3.1.2  Reference  System  Model 

At  this  point  of  the  report,  the  complete  state  vector  and  measure- 
ment equations  are  defined  only  for  the  optimum  filter/reference  system. 
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The  description  to  be  presented  here  will  clarify  the  relationships  of 
the  various  models  presented  in  Section  1. 

Table  3.  1 contains  the  reference  system  states  used  in  the  optimum 
filter.  Only  four  satellites  are  contained  in  the  state  vector  at  any  point  in 
time.  This  has  been  done  to  reduce  computer  costs.  The  four  satellites 
contained  in  the  state  vector  are  those  which  are  currently  being  tracked 
by  the  receiver.  As  the  set  of  satellite  changes,  appropriate  changes 
must  be  made  to  the  reference  system  /optimum  filter  covariance  matrix. 
This  is  analytically  justified  since  covariance  values  associated  with  un- 
tracked satellites  do  not  affect  other  system  states  either  in  the  time 
propagation  or  in  the  update.  The  only  case  where  a problem  would  arise 
is  if  a satellite  which  is  being  tracked  is  replaced  for  a short  time  and 
then  reacquired.  Should  this  situation  arise,  the  state  vector  size  will 
have  to  be  increased  unless  it  is  determined  that  ignoring  the  correlation 
introduces  negligible  error. 

The  reference  system  F matrix  and  G matrix  [see  Eq.  (3.  1)]  are 
a sparse  matrices.  Because  of  this,  displaying  only  the  non-zero  ele- 
ments is  more  enlightening.  The  non- zero  elements  of  F are  given  in 
Table  3.  2.  In  Table  3.  2,  the  parameters  in  the  value  column  are  taken 
from  Section  1.  The  values  are  listed  below  for  the  typical  models  given 
in  Section  1.  The  values  in  Table  3.  2,  are  parameterized  for  ease  in 
changing  when  new  information  is  available.  The  definitions  and  values 
are  as  follows: 

a = 

to  = w (a) 1/2 

a 1 

/>  - „3  • 10+M«° 

M = gain  at  to>2  in  dB  (for  crystal  clocks  * -220  [24]) 
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Table  3.  1 


Reference  System  State  Vector 


No. 

Symbol 

Definition 

1 

Ax 

position  error  component  x in  ECI 

2 

Ay 

position  error  component  y in  ECI 

3 

Az 

position  error  component  z in  ECI 

4 

Ax 

velocity  error  component  x in  ECI 

5 

Ay 

velocity  error  component  y in  ECI 

6 

Az 

velocity  error  component  z in  ECI 

7 

Ax 

acceleration  error  component  x in  ECI 

8 

Ay 

acceleration  error  component  y in  ECI 

9 

A'z 

acceleration  error  component  z in  ECI 

10 

C1  ) 

11 

C2  / 

12 

C3  ( 

Noise  model  for  user  clock 

13 

C4  ) 

14 

C5 

frequency  offset  for  user  clock 

15 

C6 

aging  coefficient  for  user  clock  time  error 

16 

C7 

Time  error 

17 

lTl 

Tropospheric  delay  uncertainty 

18 

position  error  component  x in  ECI  of  tracked  satellite 

19 

i 

z 

position  error  component  - in  ECI  of  tracked  satellite 

20 

position  error  component  z in  ECI  of  tracked  satellite 

21 

T.J 

Time  error  of  tracked  satellite  1 

22 

Time  rate  error  of  tracked  satellite  1 

23 

Ionospheric  residual  error  along  LOS  to  satellite  1 

24-29 

same  as  18-23  for  tracked  satellite  2 

30-35 

same  as  18-23  for  tracked  satellite  3 

36-41 

same  as  18-23  for  tracked  satellite  4 

42 

XEM1 

EM-log  random  walk  error 

43 

XEM2 

EM-log  G- Markov  error 
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Table  3.  1 (cont'd) 
Reference  System  State  Vector 


No. 

Symbol 

Definition 

44 

XGCX 

GC  random  walk  error 

45 

xgc2 

GC  G-Markov  Error 

46 

*oi 

OMEGA  bias  for  received  station  1 

47 

*02 

OMEGA  G-Markov  error  for  received  station  1 

48 

X 1 1 

49 

03  1 

x 1 ( 

04  ' 

OMEGA  periodic  error  for  received  station  1 

50-53 

same  as  46-49  for  received  OMEGA  2 

54-57 

same  as  46-49  for  received  OMEGA  3 

58-61 

same  as  46-49  for  received  OMEGA  4 
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Table  3.  2 

Reference  Syetem  F-matrix 


f 


ft 


' 


Row 

Col 

Value* 

1 

4 

1. 

1 

2 

5 

1. 

3 

6 

1. 

4 

7 

1. 

5 

8 

1. 

» 

6 

9 

1. 

10 

10 

-«a 

11 

10 

a(a-  l)wa 

I 

11 

11 

Vfclj 

12 

10 

o?(a-  l)wa 

12 

11 

a3(a-  i)w 
a 

12 

12 

-•‘“a 

1 

13 

10 

(Wj  - «0)/a3 

13 

11 

(Wj  - w0)/of2 

13 

12 

(Wj  - <u0)/a 

1 

13 

13 

-w0 

14 

15 

1. 

16 

10 

/J/a3 

16 

11 

fi/at2 

1 

16 

12 

fi/a 

. 16 

13 

16 

14 

l. 

1 

21 

22 

l. 

23 

23 

-1.  /T{ 

27 

28 

1. 

29 

29 

-1.  /T. 

1 

33 

34 

1. 

*See  text  for  variable  definitione. 


Units 


-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 

-1 

sec 


-1 

sec 


-1 

sec 
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Table  3.  2 (cont'd) 
Reference  System  F-matrix 


I 

* 


Row 

Col 

Value 

Units 

35 

35 

-l./Tj 

-1 

sec 

39 

40 

1. 

41 

41 

-1.  /T. 

-1 

sec 

43 

43 

i 

-1 

sec 
sec*  * 

45 

45 

^tgc 

47 

47 

_1*  /tom 

-1 

sec 

48 

49 

1. 

49 

48 

2 

OM 

-2 

sec 

51 

51 

•UTOM 

• 1 

sec 

52 

52 

1. 

53 

52 

*“OM 

-2 

sec 

-1 

sec 

55 

55 

•Utom 

56 

57 

1. 

57 

56 

'"OM 

-2 

sec 

59 

59 

/tom 

- 1 

sec 

60 

61 

1. 

61 

* 

60 

-CO  2 
OM 

-2 

sec 

See  text  for  variable  definitions. 
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Wj  = 6«n  2/(ir  • t2) 

u>2  = ir/(2  Tjfei  2) 

f 1»  T2.  are  breakpoints  of  the  Allan  variance  curve 

T j = .5  sec  (typical  crystal  clock  [24] ) 

T 2 - 80  sec  (typical  crystal  clock  [24]) 

T ^ = 50000  sec  (typical  crystal  clock  [24]) 

T.  = 1.  17  hr. 

1 


tem 

. 15  hr. 

T s 

GC 

23  min. 

r = 

1 hr. 

OM 

60  = 
OM 

(2lf  / 12 ) hr 

In  the  reference  system  state  vector,  states  42  through  61  will 
be  used  only  when  the  Gyrocompass,  EM-log,  and  OMEGA  are  being 
measured.  In  the  bulk  of  the  NAVSTAR/GPS  analysis,  only  the  first 
forty-one  states  will  be  propagated.  This  again  is  done  to  conserve 
computer  costs. 

To  complete  specification  of  equation  (3.  1),  the  matrix  G must 
be  defined.  However,  certain  of  the  models  are  specified  only  in  terms 


71 


of  covariance  propagation.  For  these  models  it  is  easier  to  define  the 
T T 

product  G G . The  non- zero  elements  of  G G are  listed  in  Table  3.  3. 

To  get  a particular  formulation  of  Eq.  (3.  1),  any  of  the  non-unique  square 
T 

roots  of  G G may  be  used.  The  easiest  determined  square  root  is  prob- 
ably the  one  obtained  using  a Cholesky  decomposition  [20].  Some  of  the 
values  in  Table  3.4  are  in  terms  of  parameters.  Some  of  them  are  de- 
fined above,  the  remainder  are  as  follows: 

O O O , i = x,  y,  z are  PSD  values  for  random 
posi  veil  acci  7 

walk  models 


ffCLK 


is  the  PSD  for  the  satellite  clock  fractional  frequency 

~ in-11 

error  ~ 10 


0Tri  i = 1,  2,  3,  4 is  the  prediction  residual  for  the  ionospheric 

error  

<yR.  = a esc  [ /e.2  + (18°)2] 

Ea  is  the  elevation  angle  to  the  i**1  satellite 
a is  the  RMS  correction  error  (~  8 - 17  ft)[23] 


R. 

i 


i = 1,  2,  3,  4 ionospheric  pierce  points  of  i4*1  line  of 
sight  vector 


D ionospheric  distance  constant  ~2500  km 


The  model  for  the  correlated  ionospheric  error  is  discussed 
more  fully  in  Appendix 

The  final  specification  of  the  reference  system  is  the  measure- 
ment equation.  This  equation  relates  the  state  vector  elements  to  the 
measured  quantities  or  measurables.  In  terms  of  the  state  vector  of 
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Table  3.  3 

T 

Reference  System  G G 
(Upper  Triangle  of  Symmetric  Matrix) 


* 

Row  Col  Value  Units 


1 

1 

1 

a2 

2 . 

meter  /sec 

2 

2 

pos  X 

a2 

2 

meter  /sec 

k 

3 

3 

pos  y 

a2 

2 

meter  /sec 

W 

4 

4 

pos  z 

a2  , 

. 2 3 

meter  /sec 

5 

5 

vel  x 
a2  , 

2 3 

meter  /sec 

1 

6 

6 

vel  y 
a2  , 

. ?,  3 

meter  /sec 

7 

7 

vel  z 

2 

a 

. 2,  5 

meter  /sec 

8 

8 

acc  x 

a2 

. 2,  5 

meter  /sec 

1 

9 

9 

acc  y 

a2 

2/  5 
meter  /sec 

10 

10 

acc  z 

u 2(a  - l)2 

-2 

sec 

» 

10 

11 

a 

2 2 
wa  o(«  - 1) 

-2 

sec 

10 

12 

2 2 2 
w a (a  - l) 
a 

-2 

sec 

10 

13 

(a  - l)w  (w.  - wJ/a3 
a 1 0 

-2 

sec 

• 

10 

16 

(a  - Dw^fl/a3 

-1 

sec 

11 

11 

2 2 2 
w a (a  - 1) 

-2 

sec 

11 

12 

2 3 2 

(0  a (a  - l) 

-2 

sec 

> 

11 

13 

a 

C*M0l  - l)(Wj  - “J0)/0(2 

-2 

sec 

11 

16 

a»a(a  - Dp /ot2  ' 

-1 

sec 

1 

12 

12 

m 2 (a  - !)V 

-2 

sec 

> 

♦ 

See  text  for 

variable 

a 

definitions. 
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Table  3.  3 (cont'd) 

T 

Reference  System  G G 
(Upper  Triangle  of  Symmetric  Matrix) 


Row 

Col 

* 

Value 

Units 

12 

13 

uja  - 1)(oj1  - a>o)/a 

-2 

sec 

12 

16 

to  (a  - i)0/a 

a 

-1 

sec 

13 

13 

(u5l  ■ Uo)Z/ab 

-2 

sec 

13 

16 

(Wj  - wQ)0/a6 

-1 

sec 

16 

16 

fl2  / 6 

p /a 

21 

21 

a 2 

CLK 

2. 

sec  /sec 

23 

23 

"R!2  • V\ 

2. 

sec  /sec 

23 

29 

a o' 
R1  R2 

' 2 * exp( - 1 R j - R2I/D)/T. 

2. 

sec  /sec 

23 

35 

<7  a 
R1  R3 

• 2*  exp(  - |Rj  - R3l/D)/r. 

2. 

sec  /sec 

23 

41 

a a 
R1  R4 

2 exp(-|R  - R l/D)/T. 

14  l 

2. 

sec  /sec 

27 

27 

"CLK 

2. 

sec  /sec 

29 

29 

°R22  ' 2/Ti 

2. 

sec  /sec 

29 

35 

o a 

R2  R3 

2 • exp(-|R2  - R3|  /D)/T. 

2. 

sec  /sec 

29 

41 

a a 

R2  R4 

2 • exp(  - |r  - R |/D)/T. 

0 4 l 

2. 

sec  /sec 

33 

33 

"CLK2 

2 . 

sec  /sec 

35 

35 

ffR32'  2/Ti 

2. 

sec  /sec 

35 

41. 

a a 
R3  R4 

2 ' e*p{ -|r^  . R4I/D|/t. 

2. 

sec  /sec 

19 

39 

"CLK2 

2 . 

sec  /sec 

41 

41 

'R42  • 2ITl 

2. 

sec  /sec 

«f 

44 

. 03722 

fps^/sec 

9mm  Mt  tof  ••riaM*  4»linitien«. 
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Table  3.  3 (cont'd) 

T 

Reference  System  G G 
(Upper  Triangle  of  Symmetric  Matrix) 


Row 

Col 

* 

Value 

Units 

43 

43 

6. 6 X 10'4 

fps^/sec 

44 

44 

. Z82  X 10"5 

deg^/sec 

45 

45 

1. 775  X 10‘4 

deg^/sec 

47 

47 

. 00277 

centicycle^/sec 

51 

51 

. 00277 

centicycle^/sec 

55 

55 

.00277 

centicycle^/sec 

59 

59 

. 00277 

2 

centicycle  /sec 

See  text 

for  variable  definitions. 

error  states  as  presented  here,  the  measured  quantity  is  the  difference 
between  the  actual  measured  quantity  and  the  expected  measurable  based 
on  current  estimates.  This  means  that  the  measurement  equation  can 
be  linearized  about  the  current  estimate.  The  non-zero  elements  of  the 
measurement  matrix  H are  presented  in  Table  3.4.  The  measurement 
matrix  elements  described  there  are  for  the  most  general  case  being 
considered  where  there  are  thirteen  measured  quantities.  Rows  1-4 
are  concerned  with  time-of-arrival  measurements,  rows  5-8  are  con- 
cerned with  doppler  measurements,  row  9 is  concerned  with  EM-log 
measurements,  row  10  is  concerned  with  gyrocompass  measurements, 
rows  11-13  are  concerned  with  OMEGA  station  pair  measurements. 

Any  subset  of  these  measurements  may  be  used  in  a specific  scenario. 

For  most  cases  only  the  first  four  or  eight  measurables  will  be  considered. 
A more  detailed  explanation  of  the  measurement  equations  follows  the 
symbol  defintions. 

The  measurement  matrix  is  a truly  time  varying  matrix  which 
is  also  dependent  on  the  scenario.  For 'this  rj^fson  most  of  the  elements 
are  defined  in  terms  of  variables.  The  following  are  the  definitions  of 
the  variables  from  Table  3.4  and  the  equations  following  the  definitions. 


i 

u 

X 


u1  i = 1,  2,  3,  4 are  the  x,  y,  and  z components 
of  the  unit  vector  from  the  user  to  the  i tracked  satellite 


k1 

j 


1,  2,  3,  4 j = x,  y,  z 


= distance  from  user  to  the  i tracked  satellite 

= j**1  component  of  the  relative  velocity  of  the  user 

th 

with  respect  to  the  i tracked  satellite 


Table  3.  4 

Reference  System  Measurement  Matrix 


tow 

Col 

Value 

Units 

1 

1 

1,1. 
X c 

sec /meter 

2 

M. 

U (-) 

y c 

sec  /meter 

3 

i,l. 

z c 

sec/meter 

16 

1. 

# 

17 

1. 

18 

1.1. 

V? 

sec/meter 

19 

i,i. 
y c 

sec /meter 

20 

1,1, 
z c 

sec /meter 

21 

1. 

1 

23 

1. 

2 

1 

2.1. 
U (-) 

X c 

sec /meter 

2 

2 

2,1, 
y c 

sec  /meter 

2 

3 

2 1 

U (-) 

z c 

sec /meter 

2 

16 

1. 

See  text  for  variable  definitions 


Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

Value 

Units 

2 

17 

1. 

2 

24 

u2(— ) 
X c 

sec/meter 

2 

25 

2,1, 
y c 

sec  /meter 

2 

26 

2,1, 
z c 

sec  /meter 

2 

27 

1. 

2 

29 

1. 

3 

1 

3,1. 

ux(7> 

sec  /meter 

3 

2 » 

3,1, 
y c 

sec /meter 

3 

3 

3,1, 
z c 

sec  /meter 

3 

16 

1. 

3 

17  , 

1. 

3 

30 

u3(-) 

X c 

sec/meter 

3 

31 

3,1, 
y c 

sec  /meter 

3 

32 

3.1, 

z c 

sec/meter 

See  text  for  variable  definitions. 
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Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

Value 

Units 

3 

33 

1. 

3 

35 

1. 

■ 

4 

1 

u (-) 
X c 

sec  /meter 

4 

2 

4,1. 
y c 

sec  /meter 

4 

3 

4,1. 
z c 

sec  /meter 

4 

16 

4 

1. 

4 

17 

1. 

4 

36 

4,1, 

X c 

sec  /meter 

4 

37 

u (-) 

y c 

fie  c /meter 

4 

38 

u (-) 
z c 

sec /meter 

4 

39 

1. 

4 

41 

1. 

5 

1 

kl 

X 

counts  /meter 

5 

2 

kl 

y 

counts  /meter 

* 

See  text  for  variable  definitions. 
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Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

♦ 

Value 

Units 

5 

3 

k1 

z 

counts /meter 

5 

4 

”,  ’ “Tdt(c’ 

counts  /(me ter  /sec) 

5 

5 

u‘  • 

counts  /(meter  /sec) 

5 

6 

U1  ' u)  dt(i} 
z T c 

counts  /(meter  / sec) 

5 

10 

dt  wT  fi/ei  3 

counts /sec 

5 

11 

dt  «T0  /a2 

counts /sec 

5 

12 

dt  to>T/S  /a 

counts /sec 

5 

13 

dt  Uj,P 

counts /sec 

5 

14 

dt  a>T 

counts / sec 

5 

22 

-dt  WT 

counts /sec 

6 

1 

k2 

X 

counts  /meter 

6 

2 

k2 

y 

counts /meter 

6 

3 

k2 

counts /meter 

£ 

See  text  for  variable  definitions. 
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Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

V alue 

Units 

6 

4 

• tOdt  (— ) 
x c 

counts  /( meter  /sec) 

6 

5 

• CO  dt  (— ) 

y c 

counts  /(meter  /sec) 

6 

6 

2 , ,1» 
u • <0dt  (— ) 
z c 

counts  /(meter  / sec) 

i 

6 

\ 10 

dt  uTP  /a. 3 

counts /sec 

6 

11 

dt  coT0  /a  2 

counts / sec 

6 

12 

dt  uiTP/a 

counts /sec 

x 

6 

13 

dt  t*>T0 

counts /sec 

6 

14 

dt  WT 

. 

counts /sec 

6 

28 

-dt  coT 

counts /sec 

7 

1 

k3 

X 

counts /meter 

■ 

7 

2 

k3 

y 

counts  /meter 

| 

7 

3 

k3 

z 

counts  /meter 

I 

7 

4 

u3  • CO  dt  (^) 
x T c 

counts  /(meter  / sec) 

[ 

7 

5 

« 

u3  * to  dt  (-) 
y T c 

counts /(meter /sec) 

See  text  for  variable  definition. 

c 
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Table  3.  4 (cont'd) 


Reference  System  Measurement  Matrix 

* 


Row 

Col 

Value 

Units 

7 

6 

u3  * tol  dt  (— ) 
z T c 

counts /sec 

7 

10 

dt  wTj3  /a  3 

counts /sec 

7 

11 

dt  toT0/a  2 

counts /sec 

7 

12 

dt  coT0  /a 

counts /sec 

7 

13 

dt  «T0 

counts /sec 

7 

14 

dt  WT 

counts /sec 

7 

34 

-dt  uT 

counts /sec 

8 

1 

k4 

counts  /meter 

X 

8 

2 

k4 

counts  /meter 

y 

8 

3 

k4 

counts  /meter 

z 

8 

4 

u4  • GO  dt  (— ) 
x T c 

counts/(meter  /sec) 

8 

5 

u4  • GO  dt  (— ) 
x T c 

counts  /(meter  /sec) 

8 

6 

u4  • GO  dt  (— ) 
z T c 

counts /(meter /sec) 

8 

10 

dt  wT  /J/a3 

counts /sec 

*See  text  for  variable  definition. 
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Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

Value 

8 

11 

dt  WT0/V 

1 

8 

12 

dt  ojt0  /« 

1 

8 

13 

dt  a>T0 

8 

14 

dt  WT 

1 

8 

40 

-dt  a>T 

9 

4 

H 

X 

» 

9 

5 

H 

y 

1 

9 

6 

H 

z 

9 

42 

1. 

1 

9 

43 

1. 

10 

4 

H>l 

1 

10 

5 

H*/|y| 

10 

6 

m 

> 

10 

44 

1. 

See  text  for  variable  definition. 


Units 

counts /sec 


counts /sec 


counts /sec 


counts /sec 


counts / sec 


sec /meter 


sec/meter 


sec/meter 


I 
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Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

Value 

10 

45 

1. 

11 

1 

gll 

11 

2 

g 12 

11 

3 

g13 

11 

46 

1. 

11 

47 

1. 

11 

48 

t 

1. 

11 

50 

-1. 

11 

51 

-1. 

11 

52 

-1. 

12 

1 

g2 1 

12 

2 

g22 

12 

3 

g23 

12 

46 

1. 

* 

See  text  for  variable  definition. 
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Units 


meter 


-1 


meter 


meter 


-1 


meter 


meter 


meter 


-1 


J 


I 


Table  3.  4 (cont'd) 

Reference  System  Measurement  Matrix 


Row 

Col 

Value 

Units 

1 

12 

47 

1. 

12 

48 

1. 

• 

12 

54 

-1. 

12 

55 

-1. 

1 

12 

56 

-1. 

13 

1 

g31 

-1 

meters 

f 

13 

2 

8 32 

-1 

meters 

» 

13 

3 

g33 

-1 

meters 

13 

54 

1. 

1 

13 

55 

1. 

13 

56 

1. 

» 

13 

58 

-1. 

13 

59 

-1. 

1 

13 

60 

-1. 

See  text  for  variable  definition. 
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• 1 

P = range  rate  from  the  satellite  to  the  user 
dt  is  the  Doppler  integration  time 
<*>^,is  the  satellite  transmitter  frequency 
c is  the  speed  of  light 

H , H , H are  the  x,  y,  and  z components  in  ECI  of  the  vector 
x y z 

along  the  heading  of  the  user 

H+,  H+  are  the  x,  y,  and  z components  in  ECI  of  the  vector 

y z 

in  the  plane  tangent  to  the  Earth  at  the  user's  location 
and  orthogonal  to  the  heading  vector 


| v I is  the  magnitude  of  the  user  velocity  to.  r.  t.  the  Earth 

g. . i = 1,  2,  3 ; j = 1,  2,  3 are  elements  of-A  where 
ij  ** 

b<P, 


1 -1  0 

1 0 -1 

0 0 


0 
0 

1 -1 


A A 

dL  ax 


a<pB 

SL 


ax 


b<pc 
aL  a X 


ax 


aL 

a l 

a l 

dx 

ay 

a z 

ax 

ax 

ax 

ax 

dy 

a z 

with  the  partial  derivatives  evaluated  at  the  present  user  location 
and  from  [7] 


a<p.  R. 

Tl  = *9974  ^ 


f sinL  cos  L.  cos(X  - X.)  - cos  L_  sin  L. 
I iv  1 iv  l R 1 

' L (i  - v>1/2 
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a<p  Rf 

JT  - -9974  — 


coa  LR  cos  sin  (X  - X . ) 


(1  - uf,l/2 


Rq  = radius  of  the  Earth 


LR,  XR  are  latitude  and  longitude  resp.  of  the  user 


Ly  Xf  i = A,  B,  C,  D are  the  latitude  and  longitude  resp.  of 
the  OMEGA  transmitter 


U.  = sin  L_  sin  L.  + cos  L_  cos  L.  (cos  X„  - cos  X.) 

1 R i R i R i 


f is  the  frequency 


dL  dL  dL  SX  5A  dX  . J , . , . 

Sjp'  ' "Sz'  ’ ~dx  ’ dy  ’ Can  “e  comPutec*  from  the  relations 

for  an  ellipsoid  of  revolution. 


x = R 


cos  L sin  (Wt  + X ) 
R R 

e /i  r2  • 2 r i1/2 

(1-4  sm  LR) 


= R 


cos  cos  (Wt  + XR) 

e ~ ~Z  . 2 ' .1/2 

(l-<  sm  Lr) 


z = R 


(1-4  ) sin  Lr 

e " .2  . 2.  .1/2 

(1-4  sin  Lr) 


4 = 1/298.25 
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The  modeling  section  presents  the  error  model  dynamics  for 
each  of  the  reference  system  states.  The  following  will  show  how  each 
of  the  measurables  depends  on  these  system  states.  The  measurables 
are  of  five  different  types.  They  are:  (1)  TOA  (time-of-arrival),  (2)  dop- 
pler,  (3)  EM-log  derived  velocity,  (4)  gyrocompass  heading,  and  (5)  OMEGA 
phase  errors.  The  equation  for  each  measurable  type  is  presented  be- 
low with  a brief  explanation. 

The  measurement  equations  presented  are  all  linearized  about 
a nominal  or  estimated  value.  Thus  the  equations  presented  are  the  mea- 
surement equations  of  the  deviations  of  each  quantity  from  its  expected 
or  computed  value. 

(1)  TOA  Measurements.  (Rows  1-4)  The  equation  relating  the 
difference  in  the  measured  TOA  and  the  predicted  TOA  (i.  e.  , A TOA)  to 
the  error  state  variable  is  given  by 

i i 

u u 

ATOA.  = — Ax  + — ^ Ay 

ic  c 

i i 

+ — As1  + As1  + — As1  + T * + T.  + a v* 

exeyez  sli  I 

where  i indicates  the  quantities  are  related  to  the  i**1  tracked  satellite 
and  <7v*  is  the  receiver  measurement  error  (i  * 1,  2,  3,  4). 

It  can  be  seen  from  this  equation  that  the  TOA  error  is  a function  of 
both  time  and  position  error  states  of  the  user  and  the  transmitting  satel- 
lite. Additionally  it  is  a function  of  the  atmospheric  errors  along  the 
transmission  path. 

(2)  Doppler  Measurements.  (Rows  5-8)  The  equation  relating  the  dif- 
ference between  predicted  and  measured  doppler  counts  (frequency  integrated 
over  dt)  to  the  error  state  variables  car  be  derived  from  the  relation  given  in 
Section  1. 


+ — Az  + c_  + T 
c 7 T 


u 
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I 


• » 


u 

An*  = A dt  u>  - dt  Aw*  - — dt  Ap*  + a ___  v* 
T T c DOP  2 


The  doppler  count  error  comes  from  three  basic  sources,  the 
dt  error  which  is  a function  of  user  clock  parameters,  the  error 
which  is  a function  of  satellite  clock  parameters,  and  the  p error  which 
is  a function  of  user  position  and  velocity  errors.  Each  of  these  can  be 
identified  in  terms  of  the  reference  system  states. 

The  error  in  the  integration  time  is  equal  to  first  order  to  the 
user  clock  rate  error  times  the  nominal  integration  time  dt. 

tit  = d*  (^3  cl+f?C2+«‘3  + / ,c4  + c5) 

Similarly  the  transmitted  frequency  error  is  to  first  order  equal 
to  the  satellite  clock  rate  error  times  the  nominal  transmitted  frequency. 


A«‘  = 
T 


T _ W 
s2  T 


The  range  rate  error  is  given  in  Section  1 in  terms  of  user  position 
and  velocity  errors.  Using  the  notation  above 

A fi  = ^ ( k*  Ax  + k*  Ay  + k*  Az  ) - u Ax  - u Ay  - u Az 

WT  dt  \ x y 7 z / x y 7 z 

The  measurement  error  represented  by  ffQQpV2  *8  a truncation 
error  since  only  full  counts  are  measured.  The  full  equation  is 


An1  = 


d*  “t(5  C1  * 4 c2  c3  * ?c4  + c5)  • dt  V.'i 

vo(  a / 

u / 

+ — — dt  (u  Ax  + u Ay  + u Az  1 - k*  Ax  - k*  Ay 
c \x  y z lx  y 


- k*  Az  + O v* 
z DOP  2 
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(3)  EM  log  velocity.  (Row  9)  The  EM  log  velocity  is  measured  along 
the  ships  heading  axis.  The  difference  between  the  predicted  and  mea- 
sured velocity  is  dependent  upon  the  velocity  error  projected  on  the  ships 
heading  vector  and  the  EM  log  errors.  The  EM  log  error  measurement 
is  given  in  terms  of  the  reference  system  states: 


A.  = Hx  Ai  + Hy  Ay  ♦ HZA i + XEM1  + XEM2  + <Tv3 

(4)  Gyrocompass  derived  heading.  (Row  10)  The  difference  between 
the  gyrocompass  derived  heading  and  the  predicted  velocity  vector  (in  the 
absence  of  ocean  currents)  can  be  expressed  in  terms  of  the  reference 
system  error  states.  Figure  1 shows  the  relationship  between  the  pre- 
dicted velocity  unit  vector  and  the  measured  heading.  To  first  order  the 
heading  error  can  be  represented  by 


Figure  3.  1 
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(5)  OMEGA  phase  errors.  (Rows  11*13)  The  OMEGA  phase  dif- 
ference between  the  predicted  and  measured  values  are  a function  of  the 
OMEGA  phase  errors  and  the  user  location  error  states. 


Atp 


KJ 


do 

k j Ax  + JkAy  + 


b<p 

KJ  , K 
^ <tz  + X-, 
Oz  01 


J , K J 
X01  + X02  ' X02 


I 


. K J 

+ . - X-  . 

03  03 


S<PKJ  _ 5<PKJ  S<PKJ  b X 
b x b L bx  & X Sx 


<p  = <p  - (p 
KJ  J 


where 


K = A with  J = B,  C 


and 


K = C with  J = D 


The  remaining  definition  needed  for  the  filter  relations  of  (3.  6)  is 
the  R matrix.  This  is  the  measurement  noise  covariance  matrix.  The 
purpose  of  this  measurement  matrix  is  to  simulate  receiver  measure- 
ment errors  including  truncation  effects.  Table  V gives  the  non-zero 
elements  of  R.  As  is  apparent,  all  of  the  measurement  errors  are  un- 
correlated with  the  exception  of  the  OMEGA  phase  difference  measure- 
ment errors.  , As  above  in  the  measurement  matrix  definition,  it  should 
\ 

be  noted  here  that  not  all  of  the  measurement  states  and  hence  not  all 
of  the  measurement  error  covariance  terms  are  required  in  each  simula- 
tion. Typical  values  for  the  parameters  in  Table  3.  5 are: 
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C __  A = “ chip  ® 7 * 7 

TOA  4 4 10.23  X 106 


~ 24. 5 ns 


a * .29  counts  (uniformly  distributed  truncation  error) 

DOP 


aEM  = * 193  £p8  [l1 


ffGC  = . 15  degree  [1] 


aOM  = -01cycletl1 


1-10  0 


10-10 


0 0 1-1 


1 0 


0 0 


-1  1 


0 -1 
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Table  3.  5 

Reference  System  Measurement  Noise  Covariance 


Row 

Col 

Value 

Units 

1 

1 

'toa 

2 

sec 

2 

2 

<7® 

TOA 

2 

sec 

3 

3 

'"toa 

2 

sec 

4 

4 

'2toa 

2 

sec 

5 

5 

*2dop 

counts^ 

6 

6 

*2dop 

counts^ 

7 

7 

^DOP 

counts^ 

8 

8 

ff2DOP 

counts^ 

9 

9 

<y2 

EM 

. 2 
fps 

10 

10 

'2cc 

. 2 
deg 

11 

11 

2'2OM 

11 

12 

°2om 

Jl2 

11 

'OM 

12 

12 

2'om 

12 

13 

-ff2 

OM 

13 

12 

-2 

“aOM 

13 

13 

2<j>2 

OM 

See  text  for  typical  parameter  values. 


93 


3.  2 Suboptimal  Filters 


In  this  section,  several  suboptimal  filter  candidates  will  be  de- 
scribed. Only  the  general  form  of  the  filters  will  be  described  since  the 
exact  state  vector  can  be  chosen  only  after  sensitivity  analyses.  Conse- 
quently, any  conclusions  about  advantages  in  computational  burden  must  be 
reserved  until  after  an  appropriate  state  vector  is  chosen  for  the  filter. 


3.  2.  1 Kalman  Suboptimum  Filters 

The  form  of  the  Kalman  suboptimum  filter  is  exactly  the  same 
as  the  Kalman  optimum  filter,  i.  e.,  it  is  governed  by  the  equations 
given  in  Section  2.  1.  The  primary  difference  is  the  number  of  states  in 
the  state  vector.  In  addition,  some  of  the  models  for  a given  state  may 
be  different  than  the  model  for  the  same  state  in  the  optimum  filter. 

This  is  often  done  when  it  is  desirable  for  the  model  to  attempt  to  ac- 
comodate a variety  of  errors  in  a suboptimum  fashion  rather  than  con- 
centrating on  a particular  error.  This  type  of  filter  has  proven  very 
successful  in  a large  number  of  cases.  The  number  of  states  required 
and  which  state  models  to  incorporate  are  determined  by  sensitivity 
analyses  and  experience. 

3.  2.  2 Fixed/Scheduled  Gain  Suboptimum  Filters 

This  suboptimum  filter  considerably  reduces  the  number  of  op- 
erations required  at  each  measurement  time  by  eliminating  the  gain  cal- 
culation from Eq.  (3.  6).  The  amount  of  computation  in  the  extrapolation 
phase  is  also  reduced  since  a covariance  is  no  longer  needed.  The  equa- 
tions required  are: 

= *k,k-l-*k-l  + *SJ*k  ' nk^k.k-l5*-^ 


3 


1 


3 


3 
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and  defining  equations  for  K(k). 

K(k)  may  be  a single  fixed  matrix  of  gains  or  it  may  be  a set  of 
gain  matrices  from  which  the  appropriate  matrix  is  chosen  by  some  pre- 
determined criteria  or  by  some  external  signal  (e.  g.  , a maneuver  de- 
tector) . A common  method  for  selecting  the  gain  values  is  to  solve  the 
steady  state  Riccati  equation  for  the  covariance  P. 

FP  + PFT  - PHR_1HTP  + GGT  = 0 (3.  7)  ' 

and  then  use  the  P obtained  to  calculate  K 

K = PHT  [HPHT  + R]"1 

Note  to  do  this  some  nominal  geometry  must  be  selected.  The  value 
T 

of  GG  may  also  be  adjusted  to  improve  the  performance  of  the  filter. 

The  reduction  in  the  amount  of  computation  is  naturally  accom- 

/ 

panied  by  degraded  performance  and;a  possible  requirement  for  addit- 
ional states  to  meet  minimum  performance  requirements. 

3.3  Fading  Memory  Filters  f 2 1 T 

This  filter  is  a modification  to  the  Kalman  suboptimum  filter  in 
Section  3.1.  The  basic  idea  is  to  fade  out  information  based  on  past 
measurements  and  weigh  the  most  recent  measurements  more.  This 
is  an  attempt  to  overcome  effects  of  the  mismodelling  which  is  called 
a divergence  phenomenon.  This  is  said  to  occur  when  the  estimate  of 
the  state  becomes  inconsistent  with  the  error  covariance  predicted  by 
the  filter  equations. 

The  fading  memory  filter  is  a simple  modification  to  the  filter 
equations  in  Eq.  (3.  6).  The  change  is  to  Eq.  (3.  6c)  which  becomes 
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It  is  hoped  that  this  rather  minor  modification  to  the  suboptimum 
filter  in  Section  3.  2.  1 will  make  a noticeable  improvement  in  performance. 

3.4  tt  -fl  Filter  [22l 

This  is  a simple  filter  which  has  proven  successful  in  some  ap- 
plications. This  filter  would  apply  in  the  case  where  TOA  measurements 
only  are  made.  The  state  vector  would  then  consist  of  position  and  time 
estimates  and  rates  of  each.  The  basic  assumption  of  the  OL-P  filter  for 
multivariable  systems  is  that  the  measurements  and  the  dynamics  of  the 
sets  of  a variables  and  its  associated  rates  are  independent  and  do  not 
interact.  The  simplest  form  of  the  <K-P  filter  would  use  constant  values 
for  Ot  and  P.  Other  methods  for  choosing  Cl  and  P may  also  be  considered. 

Let  s represent  a generic  state  and  r its  rate.  In  the  naviga- 
tion problem  s might  represent  each  of  the  coordinates  and  the  clock 
bias  while  r represents  the  corresponding  rates.  The  &-P  filter  can 
then  be  described  by  the  following  equations: 

s+(k)  = s"(k)  +a(k)[y(k)  - s*(k)] 

r+(k)  = r‘(k)  + ftk)[(y(k)  - s‘(k))/T] 

where 

s+,  r+  represent  updated  values  and  s , r represent 
the  extrapolated  values  and  T is  the  time  between  updates 


• '(k)  = *+(k-l)  + Tr+(k-l) 

r‘(k)  = r+(k- 1) 

y(k)  is  the  state  measurement  at  time  k 

One  method  for  choosing  "optimal  Ot-P  gains"  is  to  use  the  Kalman 
filter  equations  for  the  two  state  systems  of  the  state  variable  and  its  rate 
where  only  a measurement  of  the  state  is  available  [22],  The  formulations 
are  shown  to  be  equivalent  in  Ref.  [22].  In  the  multivariable  case,  gains 
for  each  of  the  coordinate  components  of  the  state  are  computed  ignoring 
the  correlations.  Then  one  set  of  gains  will  serve  for  each  of  the  three  co- 
ordinates plus  one  set  for  the  user  bias. 

In  using  an  a -P  filter  in  the  NAVSTAR/GPS  system,  it  will  be 
required  to  generate  pseudo-measurements  since  the  measurables  are 
time  of  arrivals  and  not  positions.  Thus  in  addition  to  the  Ot-P  filter, 
an  algorithm,  such  as  a single-fix,  must  be  implemented  to  generate 
the  pseudo -measurements. 

By  their  nature,  a -P  filters  do  not  perform  well  when  the  rate 
variables  change  rapidly.  It  is  in  the  period  of  acceleration  that  this  type 
of  filter  must  be  examined  most  carefully.  Perhaps  increased  plant  noise 
or  a fading  memory  mechanization  will  make  the  filter  less  sensitive 
to  changing  velocities. 
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Section  4.  - COVARIANCE  ANALYSIS 


4. 0 Introduction 

The  purpose  of  this  section  is  to  describe  the  covariance  analy- 
sis task  of  the  NAVSTAR/GPS  Navigation  Analysis  and  Algorithm  Develop- 
ment Study.  This  description  includes  a summary  presentation  of 
covariance  analysis,  a discussion  of  the  particular  application  to  the 
candidate  algorithms  described  in  Section  3,  and  the  results  of  the 
analysis  runs  completed  to  date. 

The  objective  of  covariance  analysis  is  to  establish  the  expected 
navigation  accuracy  of  the  filtering  algorithms  selected  as  candidates 
for  use  in  the  operational  system.  The  accuracy  can  be  established 
only  with  respect  to  the  reference  model  and  the  assumed  statistics 
of  the  error  sources.  Many  of  the  factors  affecting  accuracy,  such  as 
user  to  satellite  geometry,  receiver  configuration,  satellite  selection 
algorithm,  update  rates,  etc.  , cannot  be  modelled  in  the  format  required 
for  covariance  analysis.  For  this  reason,  each  covariance  analysis 
run  has  associated  with  it  a scenario  which  consists  of  a particular 
choice  from  all  possible  combinations  of  these  factors.  Covariance 
analysis  cannot  tell  the  entire  story  of  accuracy  except  with  respect 
to  the  reference  system  and  a particular  scenario.  It  does,  however, 
give  a good  indication  of  the  expected  error  variance  when  sever, 
scenarios  are  analyzed. 

4.  1 Covariance  Analysis 

This  section  describes  in  a summary  fashion  the  problem  to  be 
solved  by  covariance  analysis  and  the  technique  of  solution.  The 
method  of  solution  described  herein  has  been  implemented  in  a general 
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» 


i 


# 
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covariance  analysis  program  CANOMIS  (Covariance  ANalysis  Of 
Multisensor  Integrated  Systems).  The  results  presented  in  subsequent 
sections  are  the  result  of  applying  CANOMIS  to  the  specific  problem  of 
GPS  navigation. 

4.  1.  1 The  Problem 


All  of  the  candidate  filters  in  Section  3 use  weighting  matrices 
or  gains  when  incorporating  new  data.  The  gains  are  computed  based 
on  such  things  as  the  assumed  estimation  error  covariance,  the  assumed 
measurement  noise  covariance,  the  assumed  input  noise  covariance, 
and  assumed  system  dynamics  and  measurement  process.  In  general, 
each  of  the  assumed  values  may  be  incorrect  at  least  in  part.  Addi- 
tionally, it  is  necessary  to  ignore  certain  known  error  sou.ce3  to 
reduce  the  computational  burden  in  actual  operational  mechanizations. 
The  purpose  of  covariance  analysis  is  to  determine  the  expected  error 
arising  from  each  of  these  sources.  For  purposes  of  discussion,  these 
error  sources  may  be  classified  into  the  following  categories: 

(1)  an  incomplete  state  vector, 

(2)  incorrect  system  matrices, 

(3)  incorrect  initial  state  statistics, 

(4)  incorrect  statistics  for  the  white  noise  processes, 

(5)  nonwhite  noise  in  the  plant  and/or  measurements. 

These  errors  in  the  modelling  can  lead  to  the  so-called 
"divergence  problem"  which,  loosely  speaking,  occurs  when  the  actual 
errors  between  the  true  states  and  the  estimated  states  become  incon- 
sistent with  the  assumed  filter  covariance.  The  greatest  concern  is 
of  course  when  the  errors  become  much  larger  than  indicated  by  the 
filter  covariance.  Two  versions  of  this  divergence  phenomenon  may 
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be  distinguished.  "True"  divergence  is  said  to  occur  when  the  covari- 
ance of  the  actual  error  becomes  unbounded  as  the  length  of  the  data 
span  increases.  "Apparent"  divergence  occurs  when  the  actual  error 
covariance  matrix  remains  bounded,  but  is  much  larger  in  some  sense 
than  the  filter  error  covariance.  In  each  case  the  estimate  of  the  true 
state  is  unreliable  so  that  the  behavior  of  the  estimator  ia  unsatisfac- 
tory. Of  less  concern,  but  a problem  nonetheless,  is  the  situation  in 
which  the  actual  estimation  error  is  substantially  less  than  indicated 
by  the  filter  covariance.  While  the  estimator  may  provide  acceptable 
estimates,  the  lack  of  knowledge  of  the  covariance  can  nhibit  actions 
based  on  the  response  of  the  estimator. 

To  discuss  the  problem  in  more  detail,  consider  the  following 
linear  systems,  each  of  which  describes  the  errors  about  seme  point. 


% 


= F x 
F-F 


G u 
F-F 


z_  = H_  x_  + v 
-F  F -F  -F 


(4.  1) 


x = F x + G u 
—s  S— 8 S— S 

(4.2) 

- z = H x + v 

—s  S— 8 — S 


where 

x is  the  filter  state  vector,  dimensioned  n 
— r F 

Fp  is  the  filter  plant  matrix 

Gp  is  the  filter  input  distribution  matrix 

u is  the  filter  plant  noise  vector  (unity  PSD  white  noise) 
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is  the 
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is  the 

filter  observation  vector 
filter  measurement  matrix 
filter  measurement  noise  vector 

reference  error  system  state  vector,  dimensioned 

reference  error  system  plant  matrix 

reference  error  system  input  distribution  ma>  'ix 

reference  error  system  plant  noise  vector 
PSD  white  noise) 

reference  error  system  observation  vector 
reference  error  system  measurement  matrix 
reference  error  system  measurement  noise  vector 


The  system  in  equation  (4.  1)  is  the  system  used  to  generate 
the  gains  to  be  used  in  recursive  filtering  of  the  data.  Section  3 
concerned  itself  with  the  problem  of  defining  the  relationships  used 
to  generate  these  gains.  Now  the  techniques  for  determining  the 
accui  icy  obtainable  with  each  of  the  gain  computation  methods  will 
be  presented.  The  results  to  be  presented  are  of  course  scenario 
dependent.  The  technique  of  covariance  analysis  determines  directly 
a statistical  accuracy  and  thus  eliminates  the  need  for  extensive 
Monte  Carlo  simulation. 

4.1.2  Solution  Technique 

The  desired  output  of  the  covariance  analysis  is  the  1-<T  error 
of  the  filter  estimates  and  the  sensitivity  of  each  estimated  state  to 


\ 

\ 

selected  reference  system  error  sources.  This  information  is 
embodied  in  the  estimate  error  covariance  matrix  P. 


P = E[(xa 


* * .T. 

*rK*  -*TP>  ] 


(4.  3) 


where 


x is  the  filter  state  vector  augmented  with  zero*  to 
make  the  subtraction  well  defined,  i.  e.  , 


— F 


— F 


and  0 is  a (n  - n ) zero  vector. 
8 F 


The  solution  technique  is  to  generate  the  matrix  P as  a func- 
tion of  time.  To  do  this,  the  equations  for  P must  be  available  and 

can  be  developed  as  follows.  First,  to  simplify  the  notation,  define 
* 

a new  vector  x 


—8 


* * 

X = X - X_  = X - 

—s  -S  — F -S 


-F 


(4.4) 


In  the  most  general  case,  the  dynamics  of  the  new  state  vector 


x is  different  from  both  x 
— s 


and  Xj,  since  the  difference  is  being 


propagated.  However,  in  all  of  the  filters  proposed  in  Section  3,  only 

kinematic  relations  are  considered  so  that  it  is  safe  to  assume  that 

* 

F in  (4.  1)  is  a subblock  of  F in  (4.  2)  and  thus  x has  dynamics 
F s ~s 

described  by  (4.  2).  More  explicitly 


x = x - x = F x 
— 8 — 8 — F 


ff1  °1  [xf1 

- - - -4-  - - - - - + G u 

s _ T A s s 

0 | 0 0 


= F x + G u 

8 8 


(4.5) 
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Note  that  G u flora  not  appear  liner  it  ia  an  artificial  rjuan- 
r r 

tity  uaed  only  in  the  filter  gain  generation  and  eatimate  extrapolation 
aaaumea  zero-mean  noiae. 

In  addition  to  the  dynamica  deacribed  by  equation  (4.  5),  there 

ia  a meaaurement  done  at  diacrete  timea.  At  the  time  of  a measure- 
* 

ment  xg  ia  replaced  by 


(4.6) 


where 

K ia  the  gain  matrix 

ia  the  vector  of  obaervationa  from  the  reference 

error  ayatem  z = H x + v 
—a  8—a  —a 


The  equation!  to  define  P are  now  available.  Uaing  equa 
tiona  (4.  3) , (4.  4),  (4.  5),  and  (4.  6)  the  following  are  obtained. 

(a)  Extrapolation 

PkSVk-l  Pk-1  Vk-1  +Qk  (4 


where 


■ A 


*(At-T)  G#gJ  <pT(At-T)  dr 


^t.k-1 


ia  the  atate  tranaition  matrix  for  (4.  2) 
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(b)  Estimation 


Pk  = 


K H 


PK 


I - 


K H 


“1  T 


(4.  8) 


K 

0 


r K-1  T 


where 


R is  the  measurement  error  covariance 
8 


R = Efv  v 1 
s s s ’ 


and  the  partitions  are  conformable  with  the  defined  vector  x 

—s 


ORINCON'8  covariance  analysis  program  CANOMIS  propagates 
these  equations  as  a special  case  of  the  more  general  problem.  This 
is  done  by  setting  logical  input  variables  (see  CANOMIS  description). 

4.1.3  Analysis  of  Results 

All  of  the  covariance  analysis  is  being  done  in  inertial  coordi- 
nates. Therefore  the  individual  axis  components  have  no  particular 
relation  to  the  navigation  coordinates.  A more  meaningful  output  is 
the  root  sum  square  (RSS)  error.  This  is  the  usual  Euclidean  Norm 
of  the  1-a  errors  along  the  component  axes.  All  of  the  position  and 
velocity  results  presented  will  be  in  terms  of  the  KSS  quantities. 

In  addition  to  results  considering  all  error  sources,  it  is  desir- 
able to  isolate  the  effects  of  certain  individual  error  sources.  In  this 


I 


* way  an  error  budget  can  be  made  which  quantifies  the  sources  of  error. 

Also  note  that  the  error  covariance  equations  (4.  7)  and  (4.  8)  are  linear, 

therefore  it  makes  sense  to  define  sensitivities  to  the  various  error 

contributions,  i.  e.  , the  variance  of  each  of  the  navigation  variables 

2 

O.  j = 1 , 2, . . . , m may  be  written  as  a sum  of  say  r input  error 
variances  times  a sensitivity  for  each. 

n 

of  = £ fff  a..  j = 1,  2, . . . , m i = 1, 2 r (4.  9) 

1 i=l  J J1 

where 

O.  is  a generic  error  source  such  as  white  noise  PSD,  initial 
condition  variances,  etc. 

a.,  is  the  sensitivity  of  the  j**1  variance  to  the  i**1  input 
^ variance. 

These  sensitivity  coefficients  are  useful  in  that  through  their 
use,  error  budgets  can  be  updated  without  extensive  simulation.  The 
method  for  determining  these  sensitivities  is  generally  to  run 
the  covariance  analysis  with  all  of  the  error  sources  set  to  zero  with 
the  exception  of  the  one  of  interest.  For  a certain  class  of  error 
source,  the  sensitivity  is  derivable  from  the  covariance  matrix  of 
the  complete  error  system  [see  Appendix  A]. 

4.1.4  Analysis  Scenario 

The  usefulness  of  covariance  analysis  lies  in  its  ability  to  pro- 
vide statistical  information  over  an  ensemble  of  errors.  To  extend  this 
philosophy  to  the  scenarios,  an  approximate  method  of  analysis  was 
used.  In  this  method,  a step  change  of  acceleration  variance  was 
introduced  into  the  reference  system.  This  was  done  to  simulate  the 
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effects  of  an  ensemble  of  acceleration  changes  and  to  determine  the 
expected  navigation  error.  The  approximation  comes  in  due  to  the 
fact  that  the  linearization  of  the  measurement  process  must  be  made 
about  the  nominal  or  zero  acceleration  trajectory.  For  short  periods 
of  time,  this  error  is  not  large  and  the  results  are  still  valid  for  com- 
parison results. 

The  values  of  the  reference  system  parameters  used  in  the 
covariance  analyses  are  listed  in  Table  4.  1.  Table  4.  2 contains  the 
initial  covariance  matrix  values  for  the  reference  system.  These 
values  were  used  in  all  of  the  covariance  analysis  computer  runs 
except  as  noted. 

4.1.5  Results  of  Covariance  Analysis 

In  this  study,  covariance  analysis  was  used  not  only  to  analyze 
the  performance  of  suboptimal  filters,  but  as  a design  tool.  The 
filter  development  became  an  evolutionary  process  with  the  covariance 
analysis  providing  the  data  for  decisions  on  the  filter  development. 

The  primary  emphasis  in  the  filter  development  was  the  unaided  ship 
receiver.  By  unaided,  it  is  meant  that  no  external  velocity  or  accelera- 
tion information  is  provided.  It  was  assumed,  however,  that  in  *11 
cases  pitch  and  roll  information  was  available  from  gyros. 

The  candidate  filters  contained  only  kinematically  related  states 
involving  position  velocity,  acceleration,  and  clock  states.  The  states 
also  had  plant  noise  added  to  adjust  the  gains.  Table  4.  3 is  a matrix 
indicating  which  states  are  contained  in  each  of  the  filters  for  which 
extensive  analysis  was  performed.  Table  4.4  gives  the  plant  noise 
1 -O  values.  The  filter  measurement  matrices  are  the  submatrices 
of  the  reference  system  measurement  matrix  corresponding  to  the 
included  states. 


Table  4.  1.  Reference  System  Parameter  Values 


Parameter 

Value 

T1 

. 5 sec 

User  Clock 

80  sec 

Allan  variance 

T2 

parameters  10 

T3 

5x10  sec 

M 

-220  db 

Ionospheric  error 
correlation  time 

r. 

1 

4212  sec 

Ionospheric  error* 
1 -a  value 

1.7x10  ^ sec 

Satellite  clock  random 
walk  driving  noise 
1-ff  value 

1 . x 1 0 * * sec/(sec 

EM  log  Markov  error 
correlation  time 

fEM 

1380  sec 

EM  log  Markov  error 
1-ff  value 

. 079  meters/sec 

Gyrocompass  Markov 
error  correlation  time 

tgc 

1 380  sec 

Gyrocompass  Markov 
error  1-9  value 

6.1x10  ^ radians 

Omega  error 
correlation  time 

tom 

3600  sec 

Omega  Markov  error 
1-ff  value 

. 0223  cycles 

Omega  sinusoid 
error  frequency 

.4 

1.45  x 10  cycles/i 

*This  value  is  for  the  single  frequency  receiver. 
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Table  4.  1.  Reference  System  Parameter  Values  (Continued) 


Parameter 


Doppler  count 
integration  time 

Satellite  signal 
frequency 

Time -of-ar  rival 
measurement  error 
1 -<J  value 

Doppler  measurement 
error  1-flr  value 

Gyrocompass  measurement 
error  1 -(rvalue 

EM  log  measurement 
error  1-ff  value 


Value 
. 2 sec 

9 

1.57542  x 10  cycles/sec 

25  x 10  ^ sec 

. 289  counts 

2.  62  x 10  ^ radians 

. 14  meters/sec 
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Table  4.2.  Reference  System  Initial  Covariance  Value  s 
(Non-zero  elements  of  upper  triangle) 


Row 

Column 

Value 

1 

1 

8 2 
1x10  meters 

2 

2 

8 2 
1 x 10  meters 

3 

3 

8 2 
1x10  meters 

4 

4 

100  (meters/sec)2 

5 

5 

100  (meters/sec)2 

6 

6 

100  (meters/sec)2 

7 

7 

1 (meter/sec2)2 

8 

8 

2 2 

1 (meter/sec  ) 

9 

9 

2 2 

1 (meter/sec6) 

10 

10 

3.  1666  x 10  ^ 

10 

11 

3.821 x 10"2 

10 

12 

3.  96  x 10'2 

10 

13 

1.61  x 10"2 

11 

11 

5.48  x 10‘2 

11 

12 

6. 15  x 10"2 

11 

13 

1. 37  x 10-2 

12 

12 

7.85  x 10"2 

12 

13 

1. 31  x 10'2 

13 

13 

6.92  x 10"2 

14 

14 

"18 

1.  x 10  (sec/sec) 

15 

15 

1.  x 10  2*  (sec)  2 

16 

16 

1.  x 10  12  (sec)2 

17 

17 

.44  x 10’16  (sec)2 

18 

18 

9.  (meters)2 

19 

19 

9.  (meters)2 

20 

20 

9.  (meters)2 
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Table  4.2.  Reference  System  Initial  Covariance  Values 

(Non-zero  elements  of  upper  triangle)  (Continued) 


Row 

Column 

Value 

21 

21 

. , 1 8 . .2 
1.  x 10  (sec) 

22 

22 

1 . x 10  24  (sec/ sec)‘ 

23 

23 

♦ 

23 

29 

* 

23 

35 

* 

23 

41 

* 

24 

24 

2 

9.  (meters) 

25 

25 

2 

9.  (meters) 

26 

26 

2 

9.  (meters) 

27 

27 

1.  x 10  *^(sec)^ 

28 

28 

-24 

1.  x 10  (sec/sec) 

29 

29 

* 

29 

35 

♦ 

29 

41 

* 

30 

30 

2 

9.  (meters) 

31 

31 

2 

9.  (meters) 

32 

32 

2 

9.  (meters) 

33 

33 

- 1 8 2 
1.  x 10  (sec) 

34 

34 

-24 

1.  x 10  (sec/sec) 

35 

35 

♦ 

35 

41 

♦ 

36 

36 

9.  (meters)2 

37 

37 

9.  (meters)^ 

38 

38 

2 

9.  (meters) 

39 

39 

1 . x 10  * ^ (sec)2 

40 

40 

1 . x 10'24  (sec/scc) 

♦Value  computed  using  satellite/user  geometry  in  formulas  given  in 
section  1. 2.  2.  3 with  € = 17  ft. 


Table  4.  2 


Reference  System  Initial  Covariance  Values 
(Non-zero  elements  of  upper  triangle)  (Continued) 


Row 

Column 

Value 

41 

41 

* 

43 

43 

-2 

1.61  x 10  (meters/sec) 

45 

45 

3.  721  x 10  5 (radians)^ 

46 

46 

. 534  (radians)^ 

47 

47 

2 

.31  (radians) 

48 

48 

2 

. 314  (radians) 

49 

49 

6. 6 x 10  ^ 

50 

50 

2 

. 534  (radians) 

51 

51 

2 

.31  (radians) 

52 

52 

2 

. 314  (radians) 

53 

53 

6.6x10  ^ 

54 

54 

2 

. 534  (radians) 

55 

55 

2 

.31  (radians) 

56 

56 

2 

. 314  (radians) 

57 

57 

6.6x10  ^ 

58 

58 

. 534  (radians )^ 

59 

59 

2 

.31  (radians) 

60 

60 

2 

. 314  (radians) 

61 

61 

6. 6 x 10  ^ 

2 


♦Value  computed  using  satellite/user  geometry  in  formulas  given  in 
section  1.2. 2.  3 with  € * 17  ft. 
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Table  4.  3.  Filter  Configuration  Matrix 


Filter  Type 

States 

11 -State 

1 0 -State 

8 -State 

7 -State 

a-  p 

x-position 

yes 

yes 

yes 

ye  s 

ye  s 

y-position 

yes 

yes 

yes 

ye  s 

yes 

z-position 

yes 

yes 

ye  8 

yes 

ye  s 

x-velocity 

yes 

yes 

yes 

ye  s 

ye  s 

y-velocity 

yes 

ye  s 

yes 

ye  s 

yes 

z-velocity 

yes 

yes 

yes 

yes 

yes 

x-acceleration 

yes 

yes 

no 

no 

no 

y-acceleration 

yes 

yes 

no 

no 

no 

z-acceleration 

yes 

ye  s 

no 

no 

no 

user  clock  bias 

yes 

yes 

yes 

yes 

yes 

use  r clock 
frequency  error 

. 

yes 

no 

yes 

no 

yes 

r ' 

# 

Table  4.  4. 

Filter  Plant  Noise 

1 — O'  Values 

/ 

Noise  Parameter 

Nominal  Value 

States  added  to 

a 

pos 

. 1/2 
1 meter/(sec) 

x,  y,  z position 

°vel 

3/2 

. 1 meter/ (sec) 

x,  y,  z velocity 

a 

acc 

5/2 

. 01  meter/(sec) 

x,  y,  z acceleration 

°uclk 

..-12  .1/2 
10  sec/(sec) 

bias  error 

°urate 

10  * *sec/(sec)*  ^/day  frequency  error 

In  the  following  subsections,  the  covariance  analysis  of  the 
filter  development  will  be  presented  for  several  candidate  filters. 
Following  this,  sensitivity  of  filter  performance  to  key  system  param- 
eters will  be  presented. 

4.  1 . 5.  1 Eleven-State  Suboptimal  Filter 

The  eleven-state  filter  includes  all  of  the  candidate  states. 

The  first  attempt  at  the  eleven-state  filter  was  a Kalman  formulation 
with  no  modifications  made  to  accommodate  unmodelled  disturbances 

>Jc 

in  acceleration.  The  results  of  this  are  shown  in  Figures  4.  1 and  4.  2. 

In  Figure  4.  1,  no  measurements  of  doppler  shift  are  made.  Figure  4.  2 

shows  the  results  when  doppler  measurements  are  made.  The  scenario 

in  both  cases  includes  time  between  measurements  of  one  second, 

2 

a 5 meter/sec  step  change  of  acceleration  at  20  seconds,  and  four 
satellites  used  at  each  measurement  time.  Figures  4.  3 and  4.4  show 
essentially  the  same  cases  except  that  at  each  measurement  only  one 
satellite  signal  was  used  in  a round -robin  fashion  using  the  same  four 
satellites.  Figure  4.  5 shows  the  results  with  a fading  memory  filter 
with  a boxcar  of  acceleration  uncertainty  of  5 meters/ sec^  for 
20  i t s 35.  The  fading  memory  filter  uses  an  exponential  fade  factor. 
The  exponent  is  .2  times  At  (where  At  is  the  integration  step  size) 
for  0 * t * 25  and  40  < t < 50  and  it  is  .6  for  25  < t * 40. 

Examination  of  Figure  4.  1 and  4.  2 shows  that  in  the  absence  of 
doppler  measurements,  the  filter  does  not  provide  good  velocity  infor- 
mation. The  position  information  is  not  substantially  worse  except 

*NOTE:  Each  figure  contains  time  plots.  Plot  (a)  shows  the  RSS 
position  error  standard  deviation  versus  time  and  plot  (b)  shows  the 
RSS  velocity  error  standard  deviation  versus  time.  Care  should  be 
taken  when  comparing  plots  since  not  all  of  the  scales  are  the  same. 
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where  the  acceleration  uncertainty  is  large.  This  would  be  the  case 
of  the  unaided  system.  The  application  of  a fading  memory  does  appear  to 
be  of  value  in  smoothing  the  effects  of  an  unmodelled  acceleration. 

Comparing  the  results  of  Figure  4.  3 and  4.  4 to  Figures  4.  1 
and  4.2,  respectively,  shows  that  there  is  a considerable  loss  of 
accuracy  when  the  measurements  are  taken  successively  rather  than 
in  batch.  This  indicates  that  the  four-channel  receiver  has  a consider- 
able advantage  over  a single -channel  receiver. 

4.  1 . 5.  2 Ten -State  Suboptimal  Filter 

The  ten- state  filter  is  identical  to  the  eleven- state  except  that 
the  user  clock  frequency  term  is  not  included.  The  user  clock  fre- 
quency error  is  a small  number  and  its  effect  may  be  at  least  partially 
absorbed  in  a clock  bias  which  is  modelled  as  a random  walk.  The 
introduction  of  a fading  memory  filter  also  reduces  the  error  intro- 
duced by  this  simplification. 

Figure  4.6  shows  the  results  of  the  ten-state  without  doppler 

2 

measurements.  The  scenario  is  a 5 meter/sec  boxcar  change  of 
acceleration  at  20  seconds  with  a return  to  zero  at  40  seconds. 

Figure  4.  7 is  the  same  reference  system  scenario  with  a fading 
memory  filter.  The  fade  factor  is  exponential  1 times  At.  At 
t = 25  seconds,  the  exponent  was  changed  to  2 times  At. 

The  results  of  ten- state  fading  memory  filters  with  doppler 
measurements  are  shown  in  Figures  4.  8 and  4.  9.  The  run  of 
Figure  4.  8 has  a fade  factor  with  . 5 times  At  in  the  exponent  changing 
to  1 times  At  at  t = 25  seconds.  The  run  corresponding  to  Figure  4.  9 
has  the  step  change  of  acceleration  decreasing  at  t = 35  seconds  and 
the  fade  factor  exponent  equal  to  . 2 times  At  for  0 S t s 25, 

. 6 times  At  for  25  < t * 40  and  back  to  . 2 times  At  for  40  < t £ 50. 


Figure  4.6  Ten-state  filter  without  doppl 


Figure  4.7  Ten-state  filter  without  doppler  with  fading  memory. 


Figure  4.  8 Ten-state  filter  with  doppler  and  large  fading  memory  factor. 


The  sensitivity  of  the  results  to  the  measurement  errors  is 
shown  in  Figures  4.  10  and  4.  11.  The  runs  are  the  same  as  the  one 
for  Figure  4.  9,  except  for  the  measurement  errors.  For  Figure  4.  10, 
the  doppler  error  was  reduced  to  a one-half  cycle  truncation  error 
^DOP  = • *45  cycles).  The  results  in  Figure  4.  1 1 are  for  a TOA 
measurement  error  O of  5 ns  and  a doppler  truncation  error  of 
about  22  degrees  (®pQp  = • 01732  cycles). 

Comparison  of  Figure  4.  6 and  4.  7 shows  the  advantage  to  be 
gained  when  using  a fading  memory  formulation  on  the  ten- state  filter. 
The  fade  factor  may,  however,  not  be  the  best.  The  effects  of  different 
fade  factors  for  the  case  where  doppler  measurements  are  made  can 
be  seen  in  Figures  4.  8 and  4.  9.  The  filter  with  the  larger  fade  factor 
in  Figure  4.  8 has  a better  response  to  transients  while  the  filter  with 
the  smaller  fade  factor  has  better  steady  state  properties. 

The  effects  of  the  receiver  measurement  errors  can  be  deter- 
mined by  comparing  Figures  4.  9,  4.  10  and  4.  11.  The  improvement 
in  velocity  accuracy  is  shown  by  the  change  from  Figure  4.  9 to  4.  10 
and  to  4.  11.  There  does  not,  however,  appear  to  be  any  improvement 
in  position  error  when  the  TOA  measurement  error  is  reduced.  This 
indicates  that  the  position  error  is  more  dependent  on  nonreceiver 
related  errors. 

4.  1.5.3  Eight-State  Filter 

The  eight-state  filter  models  position,  velocity,  and  two-clock 
states  only.  The  acceleration  terms  are  ignored.  It  would  be  expected 
that  the  filter  would  not  perform  too  well  if  there  are  accelerations. 

To  partially  compensate  for  this,  the  fading  memory  has  been  used 
with  different  fade  factors. 
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Figure  4.  10  Ten-state  filter  with  doppler  and  reduced  doppler  measurement  e 
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Figure  4.  1 1 ; Ten-  state  filter  with  doppler  and  minimum  measurement 


Figures  4.  12  and  4.  13  show  the  results  when  the  filter  memory 
was  not  faded.  In  these  runs,  a step  change  in  acceleration  uncertainty 
of  5 meters/sec*'  was  added  at  t = 20  seconds.  Figure  4.  12  is  for 
the  filter  with  doppler  measurements.  Figure  4.  1 3 for  the  filter  with 
TOA  measurements  only. 

The  effects  of  an  adaptive  fading  memory  are  shown  in  Figures 
4.  14  and  4.  15.  Figure  4.  14  is  the  result  when  the  exponential  fade 
factor  has  an  exponent  of  1 times  At  which  changes  to  2 times  At  at 
t = 25  seconds.  The  results  of  Figure  4.  15  are  for  a fade  factor  expo- 
nent of  . 2 times  At  for  0 s t s 25  and  40  < t s 50,  with  the  exponent 
changing  to  . 6 times  At  for  25  < t s 40. 

The  results  here  show  that  the  eight-state  filter  is  not  satis- 
factory in  the  unaided  case  unless  it  has  a fading  memory  mechaniza- 
tion. 

4.  1.5.  4 Seven -State  Filter 

The  seven- state  filter  models  position,  velocity,  and  user 
clock  bias.  This  filter  has  the  same  relationship  to  the  eight-state  as 
the  ten- state  has  to  the  eleven- state  filter.  A fading  memory  or 
some  other  compensation  will  be  required  in  absence  of  the  additional 
state. 

Results  of  seven-state  filter  runs  are  shown  in  Figures  4.  16 
and  4.  17.  For  Figure  4.  16  the  filter  incorporated  doppler  measure- 
ment. The  acceleration  uncertainty  was  a boxcar  of  magnitude 
5 meters/sec  lasting  from  t = 20  seconds  to  t = 40  seconds.  In 
Figure  4.  17,  the  results  are  for  an  adaptive  fading  memory  filter 
with  the  exponent  of  the  fade  factor  equal  to  . 2 times  At  for  0 s t ^ 25 
and  40  < t < 50  and  equal  to  . 6 for  25  < t s 40.  The  input  acceleration 
uncertainty  was  a boxcar  again  but  lasting  from  t = 20  to  t = 35  seconds. 
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Figure  4.  12  Eight- state  filter  with  doppler. 


Figure  4.  13  Eight-state  filter  without  doppler. 
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Figure  4.  15  Eight-state  filter  with  small  fading  memory  fact 


Figure  4.  17  Seven- state  filter  with  doppler  and  fading  memory. 
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The  results  show  that  a fading  memory  mechanization  is  needed 
for  the  seven-state  filter  to  perform  in  a satisfactory  fashion.  In 
fact,  comparison  of  Figures  4.  15  and  4.  17  shows  that  the  seven- state 
fading  memory  performs  better  than  the  eight-state  after  the  transient. 
This  is  due  to  the  fact  that  the  time  bias  rate  error  (the  eighth  state) 
looks  like  an  acceleration  error.  The  additional  clock  state  then  adds 
damping  to  the  filter  when  unmodelled  accelerations  are  applied. 

4.  1.5.5  a-fi  Filter 

This  is  an  eight-state  filter  which  has  position,  velocity,  user 
clock  bias,  and  user  clock  frequency  error.  The  a -0  filter  must 
have  a preprocessor  such  as  a single  fix  algorithm  to  supply  it  with 
pseudo-measurements.  The  pseudo-measurements,  position  and  user 
clock  bias,  are  filtered  by  four  two-state  filters,  three  identical  filters 
in  position  coordinates  and  one  other  for  the  clock.  This  filter  has 
no  provisions  for  incorporating  doppler  measurements. 

Figure  4.  18  shows  the  results  of  a covariance  run  with  a 
2 

boxcar  of  5 meters/sec  acceleration  uncertainty  starting  at 
t = 20  seconds  and  ending  at  t = 40  seconds.  Figure  4.  19  shows  the 
results  of  the  same  run  except  that  a fading  memory  filter  was 
employed  where  the  exponent  of  the  fade  factor  was  . 5 times  At. 

Figure  4.  20  shows  the  results  using  an  adaptive  fading  memory.  The 
boxcar  of  acceleration  was  shortened  to  extend  from  t = 20  to  t = 35. 
The  exponent  of  the  fade  factor  was  . 3 times  At  for  0 < t * 25  and 
40  < t < 50  and  it  was  . 9 times  At  for  25  < t < 40. 

The  a-fi  filter  does  not  work  well  without  a fading  memory 
mechanization.  The  adaptive  fading  memory  seems  to  give  better 
steady  state  response  and  better  transient  recovery.  The  a-j?  filter, 
due  to  its  inability  to  incorporate  doppler  measurements,  does  not 
give  good  velocity  estimates. 
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Figure  4.  18  a - fi  filter. 


4.  1.5.6  Other  Factors 


There  are  certain  other  factors  which  influence  the  filter  per- 
formance which  are  of  interest.  These  factors  are  involved  in  any 
filter  mechanization.  Two  of  these  factors,  geometry  and  update  time, 
will  be  discussed  in  this  section.  Rather  than  carrying  out  the  analy- 
sis for  all  of  the  filters,  one  particular  filter  was  chosen  as  an  exam- 
ple. The  filter  used  was  the  ten-state  filter  with  doppler  measurements. 

4. 1.5.  6.1  Effects  of  Geometry 

Much  analytic  work  has  been  done  to  study  the  geometric  effects 
of  the  GPS  determined  position.  The  common  measure  for  the  geome- 
tric error  scale  factor  is  GDOP  (see  discussion  in  Section  2).  GDOP 
is,  however,  a static  measure  based  on  a single-fix  least  squares 
position  estimate.  GDOP  does  provide  a convenient  performance 
measure  to  use  when  selecting  a satellite  constellation  for  navigation. 
Since  GDOP  was  not  designed  as  a performance  measure  for  recursive 
filtering,  it  is  of  interest  to  use  the  results  of  covariance  analysis  to 
relate  GDOP  to  the  expected  navigation  error.  It  is  also  of  interest 
to  relate  GDOP  to  the  doppler  derived  velocity  error.  Figure  4.  21 
shows  the  three-axis  RSS  1-ff  navigation  error  for  the  ten-state  fading 
memory  filter  (exponent  of  .2  times  At)  with  doppler  and  the  parameter 
values  in  Tables  4.  1 and  4.  2.  The  user  was  varied  in  space  and  time 
to  get  a range  of  GDOP  values  between  2 and  9.  The  values  represent 
an  approximate  steady  state  error  with  no  Acceleration  uncertainty 
in  the  reference  system. 
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Figure  4.21  Geometric  effects  on  navigation  errors. 


velocity  error  standard  deviation 
(meters/ sec) 


The  results  indicate  that  the  RSS  position  error  1 -a  value  is 
roughly  proportional  to  GDOP.  The  RSS  doppler  determined  velocity 
error  1-Cf  value  however  appears  to  be  independent  of  geometric  fac- 
tors. This  is  probably  due  to  the  fact  that  the  dominant  errors  in  the 
determined  velocity  are  time  and  frequency  factors  and  not  dependent 
on  geometry. 

» 

4.  1.5.  6.  2 Measurement  Rate 

In  the  filter  mechanization,  the  question  arises  as  to  how  often 
measurements  should  be  taken.  Naturally  the  more  frequent  the 
measurements,  the  better  the  accuracy  will  be  up  to  certain  bounds. 
However,  more  frequent  measurements  require  more  processing 
of  data.  This  then  is  an  area  of  tradeoff  since  requiring  high  data 
processing  rates  implies  either  a faster  processor  is  required  or  the 
processing  of  other  functions  will  suffer.  The  sensitivity  of  filter  per- 
formance to  update  time  is  thus  of  interest. 

Figures  4.  22  and  4.  23  show  the  results  for  a state  filter  with 
update  rates  of  2 seconds  and  . 5 seconds.  This  is  the  same  filter 
as  the  one  used  to  generate  the  results  shown  in  Figure  4.  9. 

In  Figure  4.  24  the  sequential  channel  filter  is  shown  with  the 
update  rate  increased  to  four  measurements  per  second.  This  is 
roughly  equivalent  to  a four  satellite  measurement  taken  every  second. 
This  result  can  then  be  compared  to  the  results  shown  in  Figure  4.  2 
which  is  the  same  filter  with  measurements  of  four  satellites  taken 


in  batch  each  second.  The  result  also  is  to  be  compared  with  the 
results  shown  in  Figure  4.  4 which  is  the  sequential  channel  filter 
updated  once  a second  or  a total  of  four  seconds  for  an  entire  round 
robin. 


Figure  4.  22  Ten-state  filter  with  2 sec  update  rate. 
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4.  1.5.7  Error  Source  Senflitivitiea 

The  actual  values  in  an  error  sensitivity  are  dependent  upon 
many  parameters  of  the  filter  mechanization.  The  idea  in  this  sub- 
section is  to  find  the  approximate  sensitivity  to  the  major  reference 
system  error  sources.  The  sensitivity  to  measurement  error  and 
geometry  have  been  discussed  previously.  The  error  sources  to  be 
considered  are  the  ionospheric  delay  error,  the  satellite  position  errors, 
the  satellite  clock  errors,  measurement  errors,  and  the  user  clock  error 
mode'.  The  percent  of  total  error  is  given  in  Table  4.  5 for  the  steady 
state  RSS  position  and  velocity  error  variances  for  two  particular  filters. 
The  two  filters  are  the  ten-state  and  the  seven-state  fading  memory  filters 
with  a fade  constant  exponent  equal  to  . 2 times  At. 

4.1.6  Summary 

A summary  of  the  covariance  analysis  runs  showing  the  exact 
standard  deviation  values  for  steady  state  (just  prior  to  the  unmodelled 
acceleration)  and  the  exact  standard  deviation  values  of  the  peak  transient 
is  presented  in  Table  4.  6. 

4.  2 Monte  Carlo  Verification 


The  covariance  analysis  method  employed  is  only  approximate 
since  it  attempts  to  simulate  an  entire  ensemble  of  maneuvers  while 
using  the  linearized  model  which  assumes  no  maneuver.  The  error 
made  in  this  approach  was  assumed  a priori  to  be  small.  It  is  neces- 
sary then  to  verify  this  assumption.  This  will  be  done  by  Monte  Carlo 
simulation  runs  over  selected  scenarios.  It  will  not  be,  however,  an 
extensive  Monte  Carlo  analysis. 

4.  2.  1 Monte  Carlo  Simulation 

The  Monte  Carlo  simulation  will  be  done  with  the  same  com- 
puter analysis  program,  CANOMIS,  that  was  used  in  the  covariance 
analysis.  This  ensures  that  the  same  models  are  being  used.  The 
difference  is  that  instead  of  propagating  a covariance  matrix,  a 

145 


MM 


SSMMHM 


Table  4.5  Filter  Error  Budget 


10-State  Filter 

7 -State 

Filter 

Error  Source 

% of  3-axiB 
position  error 
variance  sum 

% of  3-axis 
velocity  error 
variance  sum 

% of  3-axis 
position  error 
variance  sum 

% of  3 -axis 
velocity  error 
variance  sum 

Ionospheric  delay 
error 

54.7 

~ 0 

53.  5 

~ 0 

Satellite  position 
error  (all  satellite) 

25.  7 

~ 0 

25.  2 

~ 0 

TOA  measure- 
ment error 

19.  2 

~ 0 

18.  8 

~ 0 

Acceleration 

disturbance 

~ 0 

~ 0 

~ 0 

47. 

User  clock 

~ 0 

1. 

2. 

44.  7 

Doppler  measure- 
ment error 

- 0 

98. 

~ 0 

8. 

Table  4.6  Covariance  Analysis  Summary 


Steady  State  Standard  Deviation  Peak  Standard  Deviation 


Figure  No. 

(Extrapolation  Value) 

during  Transient 

Position 

(Meter) 

Velocity 

(Meter/sec) 

Position 

(Meters) 

V elocity 
(Meters  /sec) 

4.  1 

16.27 

2.96 

76.29 

35.23 

4.2 

13.2 

.28 

21.40 

17.  94 

4.  3 

24.68 

5.  12 

134. 67 

45.  78 

4.4 

15.92 

.62 

50.  02 

28.  C8 

4.  5 

14.  03 

.33 

15.  77 

12.  94 

4.6 

16.27 

2.98 

110.46 

38.94 

4.  7 

18.  47 

4. 18 

37.  58 

25.26 

4.8 

16.  78 

.62 

17.  33 

8.63 

4.9 

13.  51 

.33 

15.  31 

13.  08 

4.  10 

13.48 

.20 

15.29 

10.69 

4.  11 

14.  58 

.05 

15.  21 

8.  66 

4.  12 

12.99 

.64 

58.2 

28.28 

4.  13 

15.03 

1.79 

489. 90 

116.  00 

4.  14 

17.01 

1.92 

37.45 

26.  12 

4.  15 

13.70 

.67 

20.92 

20.  51 

4.  16 

13.  32 

1. 12 

557.79 

96.  76 

4.17 

14.25 

.82 

27.  03 

26.06 

4. 18 

15.  54 

1.82 

496.  68 

115.  35 

4.19 

19.36 

3.21 

55.  86 

38.  38 

Table  4.  6 (continued) 


4.20 

17.  14 

2.  02 

64.  82 

43.24 

4.22 

15.  14 

. 46 

25.  18 

19.24 

4.  23 

13.  44 

.25 

14.  39 

9.42 

4.24 

13.94 

. 52 

20.  18 

15.41 
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reference  system  state  vector  is  propagated  and  the  measurement 
matrix  is  linearized  about  a simulated  trajectory. 

The  Monte  Carlo  simulations  use  gains  generated  by  the  can- 
didate filter  and  with  these  gains  process  the  simulated  data  according 
to  equation  (4.  6).  The  simulated  data  is  generated  from  equation  (4.  2) 
by  replacing  the  noise  vectors  with  white  Gaussian  pseudo- random 
number  vectors  with  the  appropriate  variances. 

The  scenarios  used  are  the  ones  outlined  in  Section  1.  The 
satellite  selection  was  done  using  the  approximate  volume  minimiza- 
tion algorithm  described  in  Section  2.  The  verification  runs  were 
made  with  the  ten-state  filter,  the  eight-state  filter,  and  the  Of -p 
filter,  all  using  doppler  measurements  but  otherwise  unaided.  These 
were  chosen  as  representative  of  all  the  filters  analyzed. 

As  in  Section  4.  1,  the  results  will  be  shown  via  plots  of  the 
3-axis  RSS  position  and  velocity  error.  The  (a)  part  of  each  figure 
is  the  position  error  and  the  (b)  part  is  the  velocity  error. 

4.  2.  2 Monte  Carlo  Results 


The  results  of  several  Monte  Carlo  simulations  are  shown  in 
Figures  4.  25  through  4.  31.  Figures  4.  25  and  4.  26  show  the  results 
of  the  eight- state  fading  memory  filter  with  the  fade  constant  exponent 
equal  to  . 33  times  At.  Figure  4.  25  shows  the  results  for  the  ship 
scenario  and  Figure  4.  26  for  the  aircraft  scenario.  The  results  for 
a ten- state  fading  memory  filter  with  the  fade  constant  exponent 
equal  to  . 33  times  At  applied  to  the  ship  scenario  are  shown  in 
Figure  4.  27.  Figures  4.  28  to  4.  30  are  a series  of  ten- state  fading 
memory  filters  with  different  fade  constants  (exponents  equal  to 
. 2 At,  . 33  At,  and  . 5 At,  respectively)  applied  to  the  aircraft 


149 


(small  fade  factor). 


(medium  fade  factor). 


Figure  4.  30  Ten-state  fading  memory  filter  with  aircraft  scenario 
(large  fade  factor). 


scenario.  The  results  may  be  validly  compared  since  the  pseudo- 
random number  generator  was  started  from  the  same  value  in  each 
case.  Finally,  Figure  4.  31  shows  the  results  using  a fading  memory 
oc-fi  filter  with  the  fade  constant  exponent  equal  to  . 33  times  At. 

The  results  from  the  two  eight-state  filter  runs  (Figures  4.  25 
and  4.  26)  confirm  the  covariance  analysis  results.  The  mean  RSS 
position  error  in  Figure  4.  25  (a)  appears  to  be  below  the  value  pre- 
dicted by  the  covariance  analysis.  The  velocity  error  shown  in 
Figure  4.  25  (b)  is  within  the  limits  predicted  by  the  covariance  analy- 
sis. The  one  notable  deviation  in  the  velocity  error  is  at  t = 260, 
where  the  ship  undergoes  a 3 degree/sec  turn  at  15  knots.  This  turn 
is  an  acceleration  of  approximately  . 04  g so  that  the  velocity  error 
is  to  be  expected.  The  velocity  errors  shown  in  Figure  4.26  (b)  are 
for  the  aircraft  scenario.  The  aircraft  turns  at  t = 300  and  t = 360 
are  approximately  5 g and  . 5 g turns,  respectively.  The  velocity 
errors  in  these  turns  are  again  what  the  covariance  analysis  would 
predict  when  the  acceleration  uncertainties  are  scaled  to  the  values 
in  the  scenario.  The  RSS  position  error  from  Figure  4.  26  (a)  is 
approximately  at  the  1-0  value  from  the  covariance  analysis,  except 
during  the  5 g maneuver.  The  unusual  curve  shape  starting  just  prior 
to  t = 500  is  due  to  the  fact  that  the  acceleration  goes  to  zero  there. 

The  results  for  the  ten-state  fading  memory  filter  applied  to 
the  ship  scenario  (Figure  4.  27)  show  good  position  results.  The 
velocity  estimation  error  is  somewhat  larger  than  the  steady  state 
1-0  values  predicted  by  the  covariance  analysis.  The  results  of  the 
aircraft  case  (Figures  4.  28  through  4.  30)  give  results  which  are  again 
in  good  agreement  with  the  covariance  analysis.  Variation  of  the  fade 
factor  shows  the  larger  fade  factor  reduces  the  magnitude  of  the  tran- 
sient errors  during  periods  of  large  acceleration  while  the  steady  state 


errors  increase.  This  again  is  the  result  which  would  be  expected 
after  examining  the  covariance  analysis  results.  The  shape  of  the 
velocity  error  curves  indicate  the  ten- state  filter  is  able  to  track 
accelerations  after  an  initial  transient.  This  can  be  seen  from  the 
small  "bumps"  on  the  curve  at  t = 360  and  t = 480. 

The  last  Monte  Carlo  results  shown  (Figure  4.  31)  confirm 
good  agreement  between  the  covariance  analysis  and  the  simulation 
for  the  fading  memory  Oc-fi  filter. 

In  the  simulation  runs,  the  fading  memory  filters  were  mechanized 
with  the  same  fade  factor  throughout.  An  analysis  into  making  the  fade 
factor  a function  of  the  measurement  residuals  in  probably  worth  while 
as  a future  task. 


Section  5.  - COMPUTER  PROGRAM  SIZING 


5. 0 Introduction 

This  section  addresses  the  final  task  of  the  algorithm  develop- 
ment study.  This  task  was  to  determine  the  computational  require- 
ments (i.  e.  , arithmetic  operations  and  storage  locations)  for  the 
candidate  algorithms.  These  results  will  provide  an  input  to  the  effort 
establishing  the  computer  program  size.  The  results  to  be  presented 
here  do  not  account  for  computer  word  length.  The  true  computer 
size  requirement  will  probably  depend  on  the  word  length  used  in  the 
computer.  This  problem  has  been  considered  to  be  beyond  the  scope 
of  this  study,  though  ORINCON  is  aware  that  it  must  be  addressed 
somewhere. 

The  actual  computer  mechanization  of  the  navigation  equations 
can  be  divided  in  two  parts.  The  first  part  is  the  preprocessing  of  the 
raw  measurement  data  to  obtain  the  measurables  for  the  filter.  The 
second  part  is  the  filtering  of  these  measurables  to  obtain  the  actual 
navigation  information.  It  is  the  latter  problem  which  will  receive  the 
most  attention  in  this  section.  The  first  part  is  more  dependent  on 
factors  outside  this  study. 

5.  1 Common  Requirements 

All  of  the  candidate  filters  fall  into  the  category  of  extended 
filters.  By  this  it  is  meant  that  the  filter  is  linearized  about  the 
current  estimate.  The  measurables  for  this  type  of  filter  are  the 
difference  between  the  predicted  value  of  the  measurement  and  the 
actual  measured  value.  The  computation  of  the  expected  measure- 
ment is  a common  requirement  of  all  of  the  filter  mechanizations. 


The  following  steps  are  required  to  perform  the  measurable  computa- 
tion: 

(a)  Determination  of  the  time  of  transmission  of  the 
received  signal.  This  can  be  accomplished  from 
knowledge  of  the  system  requirements  which  specify 
the  times  of  transmission. 

(b)  Determination  of  the  satellite  positions  at  the  time 
of  signal  transmission.  This  can  be  done  using  the 
equations  for  this  purpose  presented  in  Section  1. 

(c)  Correction  of  the  measured  time  for  satellite  clock 
errors  and  ionospheric  delay  errors.  This  is  done 
by  using  the  received  data  and  the  algorithms  for 
these  purposes  presented  in  Section  1. 

(d)  Computation  of  expected  measurement  values.  This 
computation  involves  using  the  estimate  of  the  user 
position  and  velocity  at  the  time  of  reception  and  the 
computed  satellite  locations.  From  these,  expected 
range  and  range -rate  values  can  be  computed. 

(e)  Formation  of  the  filter  measurables.  This  final 
step  is  done  by  differencing  the  expected  measure- 
ments and  the  actual  corrected  measurements. 

Other  common  requirements  include  the  alert  computation  "> 

(such  as  the  candidates  outlined  in  Section  2),  general  matrix  operation 
subroutines,  and  special  mathematical  functions  (e.g.,  square  root, 
vector  norm,  etc.  ). 
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5.  2 Filter  Requirements 

In  this  subsection,  the  computer  requirements  for  a general 
Kalman  filter  and  the  specific  candidate  filters  will  be  presented. 

This  includes  the  storage  requirements  and  the  operation  counts. 

5.  2.  1 General  Kalman  Filter  [25] 

The  computational  requirements  for  the  Kalman  filter  are 
given  in  Table  5.  1.  Here  n is  the  number  of  states  and  m is  the 
number  of  measurables.  This  does  not  include  the  computation  for 
guaranteeing  symmetry  (Item  15).  The  computation  time  is  repre- 
sented by  the  number  of  multiplications  since  it  is  largely  governed 
by  the  number  of  multiplications.  Additions  and  subtractions  normally 
require  far  less  time  in  most  computers.  In  any  case,  these  opera- 
tions are  of  the  same  order  as  the  number  of  multiplications  so,  for 
purposes  of  comparing  algorithms,  this  operation  count  is  considered 
to  be  valid. 

It  will  be  noted  below  that  for  the  simple  dynamics  assumed 

in  each  of  the  candidate  filters  that  items  4 and  6 may  be  done  more 

efficiently.  Implementation  of  a fading  memory  filter  adds  up  to  an 
2 

addition  n multiplications  to  item  6. 

5.2.2  Eleven-State  Suboptimum  Filter 

For  the  eleven -state  filter  the  simple  dynamics  mean  that 
steps  4 and  (>  in  Table  5.  1 each  require  10n  multiplications  instead 
of  n^.  Also  there  is  no  requirement  for  storage  for  the  state  transi- 
tion matrix.  The  measurement  noise  covariance  matrix  is  also 
diagonal  so  that  only  m storage  locations  are  needed.  Another 
consequence  of  the  uncorrelated  measurement  errors  will  be  dis- 
cussed below  in  Section  5.  3. 
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Table  5.  1.  Computational  Requirements  for  the 
Kalman  Filter 


Quantity 

Multiplications 

Storage 

1. 

A 

n 

2. 

PK-1 

2 

n 

3. 

*k,k-l 

2 

n 

4. 

^k,k-lPk-l 

3 

n 

Store  in  P,  , 
k- 1 

5. 

°k 

2 

n 

6. 

pk"fk.k.ipk-1®kT.k.i+Qk 

3 

n 

Store  in  $P 
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Hk 

mn 

8. 

pxT 

2 

mn 

Store  in  P, 
k 

9. 

Rk 

2 

m 

10. 

Hkpk^+Rk 

2 

m n 

2 

m 

11. 

,HkpkHkT+Rk>'' 

3 

m 

T 

Store  in  HPH  + R 

12. 

Kk  • pkH^(Hkp^+Rk>'1 

2 

m n 

mn 
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^k 

m 

14. 

il 

** 
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mn 

Store  in 

15. 

Pk/k  = Pk/k-rKk(Pk/k-lHk) 

2 

mn 

Store  in  P 

k 

16. 

Scratch  storage 

2 

n 

TOTALS 

2n  +?.mn 

3A,  2 . 

+m  +2m  n+mn 

2 2 
4n  +n+2mn+2mn  +m 
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5.2.3  Ten-State  Suboptimum  Filter 

In  the  ten-state  filter,  steps  4 and  6 of  Table  5.  1 reduce  to 
9n  multiplications  each  with  again  no  requirement  for  a state  transi- 
tion matrix.  The  computational  requirements  for  this  filter  are  in 
Table  5.  2. 

5.2.4  Eight-State  Suboptimal  Filter 

The  simple  dynamics  used  in  the  eight-state  filter  mean  that 
steps  4 and  6 of  Table  5.  1 reduce  to  4n  multiplications.  The  results 
are  presented  in  Table  5.  2. 

5.  2.  5 Seven-State  Suboptimal  Filter 

The  transition  steps  for  the  seven-state  filter  given  in  steps 
4 and  6 of  Table  5.  1 require  only  3n  multiplications  and  no  state 
transition  matrix  storage.  The  results  are  given  in  Table  5.  2. 

5.  2.  6 Ot-0  Filter 

The  a-/9  is  an  eight- state  filter  which  is  composed  of  four 
two-state  filters.  The  transition  steps  reduce  to  8 multiplications 
and  there  is  only  one  observable  per  filter.  The  results  given  in 
Table  5.  2 do  not  include  the  prefiltering  of  the  data. 

4.  1 Additional  Considerations 

In  this  subsection,  additional  considerations  which  impact 

■i  - "■  s<.ii*'»mI  requirements  are  discussed. 


I*  • 


V.  " 

v c 

~ t m 

0-  — c 

G.  § 

0 m 

Q ” " 

^ c ^ 

T3  0)  « 

C C £ 
<i  c a 


M’ 

cm 

"t 

00  M4 

o 

00 

NO 

00 

in 

oo 

00  M* 

O' 

in  » 

vO 

o 

no  co 

CM 

o * 

cm 

co 

cm 

cm  «-H 

CM 

NO 

CM 

o 

*4 

00 

w— « 

in 

O' 

oo 

4 

co 

CM 

CM 

c 

CM 

c 

CM 

C 

CM 

C 

e 

4 

c 

4 

c 

4 

G 

4 

£ 

£ 

£ 

£ 

4 

c 

c 

c 

c 

C 

CM 

CM 

CM 

CM 

CM 

+ 

+ 

4 

4 

4 

£ 

£ 

£ 

£ 

CM 

CM 

CM 

CM 

CM 

4 

c 

+ 

4 

4 

+ 

4 

c 

G 

c 

c 

CM 

4 

4 

4 

4 

c 

<M 

CM 

CM 

CM 

co 

C 

C 

C 

c 

CO 

CO 

CO 

CO 

4 

» 


t 


I 


I 


» 


» 


I 


5.  3.  1 Sequential  Processing 

Since  the  measurement  errors  in  the  GPS  receiver  are  uncor- 
related (i.  e.  , the  R matrix  is  diagonal)  the  data  may  be  processed 
sequentially  even  though  it  is  received  in  batch  [26]  with  exactly  the 
same  results.  Processing  the  data  sequentially  always  requires 
fewer  multiplications  than  batch  processing  [25].  For  scalar  sequen- 
tial measurements,  the  totals  in  Table  5.  2 (except  the  Ct-P  filter  which 
already  uses  this  advantage)  can  be  reduced  by  [25] 

3 2 

AM  = m - m + 2nm  - 2nm 
= m(m- 1 )(m+2n) 

The  numbers  in  parentheses  in  Table  5.  2 are  the  number  of 
multiplications  reduced  by  AM.  Since  the  measurement  errors  are  all 
uncorrelated  in  the  GPS  environment,  the  sequential  processing  should 
be  implemented. 

5.3.2  Computational  Form  of  the  Filter 

Much  work  has  been  done  in  other  studies  on  the  computational 
form  of  the  filter.  Selection  of  the  form  may  have  considerable  impact 
on  numerical  errors.  The  most  common  numerical  error  seems  to 
be  that  the  filter  covariance  matrix  loses  its  positive  definiteness. 

One  method  of  solving  the  numerical  problems  of  the  covariance 
matrix  is  to  use  the  "stabilized"  form  of  the  filter.  This  is  done  by 
replacing  step  15  of  Table  5.  1 by  the  equivalent  form  given  in  equa- 
tion (3.  6d)  of  Section  3.  This  form  of  the  filter  requires  considerably 
more  multiplications.  Another  stabilization  procedure  is  to  forct- 
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the  result  of  the  covariance  reset  in  step  15  of  Table  5.  2 to  be  sym- 
metric by  averaging  the  off-diagonal  elements,  i.  e., 

p..  + p.. 

p =p  = — -Ji 

Pij  Pji  2 

\ ■ 


where 


p..,  p..  are  the  new  values 
U 


p..,  p..  are  the  results  of  step  15. 

»J 


> i 


Both  of  these  techniques  do  not  necessarily  solve  the  precision 
problem. 

A different  approach  is  through  the  square  root  formulation 
of  the  filter.  Carlson  [26]  and  others  have  developed  square  root 
fo.  r^tltions  which  solve  the  positive  definiteness  problem  and  also 
alleviate  precision  problems.  These  methods  are  well  documented 
in  the  literature  [26,27].  Another  covariance  factorization  technique, 
the  U-D  filter,  which  is  not  quite  so  well  known  is  presented  in  refe- 
rence 28.  This  formulation  is  equivalent  in  numerical  performance 
to  the  other  square  root  formulations,  but  it  does  not  require  a square 
root  to  be  taken  at  each  point  of  propagation.  For  a small  processor, 
this  cai?  lead  to  a significant  time  savings. 

A summary  of  the  U-D  filter  equations  for  the  case  of  scalar 
measurements  is  presented  here.  The  reader  is  referred  to  reference 
28  for  a thorough  explanation  of  the  U-D  mechanization. 

Suppose,  the  n-dimensional  error  covariance  matrix,  P , is 
factored  such  that 


P = UDU 


(5.  1) 
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where  U is  upper  triangular  with  unity  diagonal  elements  and 

D * diag(d  , . . . , d ).  The  matrices  U and  D are  referred  to  as 
i n 

the  U-D  factors  of  P . They  are  unique,  provided  that  P is  positive 
definite,  and  can  be  constructed  using  a Cholesky  factorization  [20]. 


U-D  Measurement  Update  Algorithm 

Given  a priori  covariance  factors  U and  D and  scalar  measure- 

2 

ment  z = Hx  + v,  where  E(v  ) = r,  the  updated  U-D  covariance  factors 
and  the  Kalman  gain  (U,  D and  K respectively)  can  be  obtained  as 
follows: 

fT  = HU  ; fT  = (f. f ) (5.2) 

i n 


v = Df  ; v.  = d.f. 

J J 1 

-T 

Kj  = (Vj , 0, . . . , 0) 


“l  ’ ' + V, 


dj  = (r/Oj)dj 


For  j = 2 n cycle  through 


a a (»  + \r  f 

J j-1  jj 


d J ■ 


XJ  ■ -VaJ-l 


(5.3) 

(5.4) 

(5.5) 

(5.6) 

(5.7) 

(5.8) 

(5.9) 
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U.  - U.  + l.K.  . 
J J J J-1 


K.  = K.  , + v.L. 
J J*1  J J 


(5. 10) 

(5.  ID 


where  U = [U_,  U_...,U  ].  The  component  U vectors  have  the  form 
12  n 


U.  = (U.(l),  — U.( j-1).  1,0,.  ..,0) 
J } J 


and  D = diag  (d, , . . . , d ).  The  Kalman  gain  is  given  by 
1 n 


K = K /a 
n n 


(5. 12) 

The  salient  feature  of  this  algorithm  is  the  way  in  which  the  updated 
diagonal  D elements  are  computed.  Since  the  quantities  are  calculated 
as  positive  sums,  it  follows  that  the  updated  d's  are  fractions  of  their 
a priori  values  and  therefore  cancellation  errors  present  in  the  conven- 
tional measurement  updating  equation  are  avoided.  The  positivity  of  D,  and 
hence  of  P,  is  therefore  assured.  Furthermore,  the  elements  of  D can 
diminish  to  near -zero  without  affecting  the  stability  of  the  algorithm. 

Modified  Givens  techniques  can  be  employed  to  accomplish  time 
updating  of  the  U-D  factors,  and  the  resulting  algorithm  is  the  following: 
Let 


[wr”2 

Wn+k] 

> 

Diag  (dj 

, . . . , ^n+j^)  J k = no.  of  columns  of  G 

— T 

3 

WDW 

can  be  computed  as  follows:  For 

j = n, ....  1 cycle  through  the  following  as  indicated, 
m:  = j + k 


(5.13) 


The  symbol 


denotes  replacement  (i.  e.  , replace  m by  j+k). 
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For  i = evaluate  Eqa.  (5-14)  - (5.27)  as  indicated. 


m m 


fit  = d.w.(j) 


d'  : = aw  (j)  + 0w  (j) 
mm  1 


» : = P Id! 


m 


d.:  = d.d  Id! 
i i m m 


(5.14) 

(5.15) 
(5.  16) 
(5.  17) 
(5.18) 


:: 

i , 


v:  = w. 
i 


(5.19) 


If  i = m-1  evaluate  Eqs.  (5.20)  - (5.23) 


c:  = a/d' 
m 


w.(4):  = w (j)v(X)  - v(j)w  (4) 
i m m 


w (4)  s = cw  (4)  + §v(4) 
m m 


w.(j):  = 0,  w (j)s  = 1 
i m 


^ = • • • i J*1 


If  l<m-l  evaluate 


wt(4)}  = v(4)  - v(j)wm(4)  | 

w (4)i  = w (4)  + sw  (4) 
in  m l 


4 s 1# • * * i jal 


When  i<m-l,  w (J)  "1 
m 


(5.20) 

(5.21) 

(5.22) 

(5.23) 


(5.24) 

(5.25) 
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w.(j)j  = 0 

i 

(5.26) 

d : = d* 
m m 

(5.27) 

Upon  completion  of  this  recursion  the  W and 

D arrays  contain  U and 

stored  as  follows 

k n 

W = [T|U][n 

(5.28) 

= V j=1 n 

(5.29) 
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SECTION  6.  - CONCLUSIONS 


6.  0 Conclusions  of  the  Study 

The  purposes  of  this  study  were  to  develop  and  analyze  navigation 
algorithms  for  use  in  Navy  NAVSTAR/GPS  navigation  receivers.  By 
performing  algorithm  design  and  analysis  for  various  receiver  configurations, 
it  was  possible  to  provide  data  useful  in  certain  trade-offs.  This  report 
has  presented  much  data  and  many  results  from  which  the  reader  may  derive 
his  own  conclusions  . Changing  requirements  and  hardware  specifications 
may  at  some  point  in  time  invalidate  certain  of  the  results  summarized 
below.  Nonetheless  the  following  conclusions  have  been  drawn  by  ORINCON. 

A.  The  Alert  algorithm  based  on  the  approximate  volume  max- 
imization or  exhaustive  search  should  be  used. 

B.  The  successive  linearizations  of  the  measurement  matrix 
method  should  be  used  to  provide  single  fixes  or  pseudo -measure- 
ments . 

C.  The  multi- channel  receiver  gives  better  transient  response 
than  the  single  channel  receiver  if  the  frequency  of  measurements 
is  the  same.  The  single  channel  performance  is  nearly  equivalent 
to  the  multi- channel  when  the  data  rate  is  equal.  This  means  that 

if  the  single  channel  receiver  can  make  n -measurements  in  the  same 
time  that  the  n -channel  receiver  makes  one  measurement,  then  the 
single  channel  performance  is  on  a par  with  the  n-channel  reciever. 
(This  conclusion  does  not  take  into  account  the  possibility  of  reducing 
ionospheric  errors  by  using  different  frequency  channels.  This  is 
.a  different  issue). 

D.  In  order  to  provide  good  velocity  estimates  from  the  GPS  data, 
doppler  measurements  must  be  made.  Without  the  doppler  measure- 
ments, position  estimates  are  only  slightly  degraded,  however  the 


velocity  errors  are  of  the  order  of  several  knots. 

E.  The  addition  of  acceleration  states  in  the  filter  provides  a 
better  velocity  estimate  in  both  the  steady  state  and  in  the 
transient  response  to  unmodelled  acceleration. 

F.  The  addition  of  a fading  memory  to  the  filter  improves  the 
performance  of  the  filter  when  unmodelled  acceleration  are  pre- 
sent. This  is  at  the  expense  of  the  steady  state  performance  when 
there  are  no  acceleration  disturbances.  Some  sort  of  adaptive 
fading  memory  would  provide  a good  compromise  solution. 

G.  The  measurement  update  rate  (within  certain  bounds)  does  not 
appear  to  be  a significant  consideration  in  systems  where  there 
are  no  unmodelled  accelerations.  The  update  rate  becomes  sign- 
ificant when  there  are  unmodelled  accelerations. 

H.  The  ionospheric  errors  appear  to  be  a significant  error  source. 
The  value  of  the  error  source  is  an  uncertain  quantity  at  this  time. 
More  data  should  be  gathered  on  the  proposed  ionospheric  cor- 
rection schemes  for  the  single  and  dual  frequency  receivers. 

I.  The  receiver  TOA  measurement  error  is  not  as  significant  an 
error  source  in  relation  to  the  overall  position  error  as  is  the 
doppler  measurement  error  in  relation  to  the  overall  velocity 
error.  This  is  of  course  only  in  the  ranges  considered  in  this 
study.  Emphasis  should  be  placed  on  improving  the  doppler  meas- 
urement accuracy. 

J.  For  low  accuracy  systems  where  an  accurate  velocity  estimate 
is  not  needed,  the  Of  -fi  filter  provides  a computationally  simple 
filter  when  combined  with  the  pseudo -measurements  generated  by 
the  successive  linearizations  algorithms. 

K.  Sequential  processing  of  the  data  should  be  done  to  save  computer 
storage  and  computation  time.  Also  a factorization  technique  should 
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I 


probably  be  used. 

I 

In  all  of  the  above  conclusions,  the  unmodelled  accelerations  refer 
to  unaided  systems.  The  amount  which  the  response  to  unmodelled  accel- 
erations should  be  weighed  in  making  any  design  trade-offs  is  of  course 
• dependent  on  the  expected  amount  of  acceleration  disturbances  (both 

frequency  and  magnitude  of  disturbances). 
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APPENDIX  A 

Filter  Sensitivity  to  a Class  of 
Unmodelled  Errors 

In  this  appendix,  a simple  method  will  be  developed  for 
attacking  a particular  problem  in  the  covariance  analysis  of  subop- 
timal  filters.  The  problem  being  considered  is  that  of  determining 
sensitivity  to  unmodeled  error  sources  of  a certain  class.  The  class 
of  unmodeled  errors  is  somewhat  restrictive,  but  it  is  rich  enough 
to  be  of  practical  importance. 

The  setting  for  covariance  analysis  is  as  follows.  First  a 
reference  system  is  defined  which  contains  all  of  the  error  sources 
modeled. 


” \k-l  -^k-1  + Fk,k-1  Hk 


(A-l ) 


whe  re 


<0, 


k,  k- 1 


is  the  n state  vector  at  time  k 

8 

is  the  n x n state  transition  matrix  from  time 
s s 
k-1  to  k 


rR  ^ j is  the  ng  x r#  input  disturbance  at  time  k 
Z, 


Xk 

^k 


is  the  m observation  vector  at  time  k 
s 

is  the  m x n measurement  matrix  at  time  k 

s s 

is  the  mg  observation  noise  vector  at  time  k 
is  the  r#  vector  of  input  noise  at  time  k 


It  is  desirable  to  estimate  certain  states  of  X^.  but  due  to 
computer  size  restrictions  certain  of  the  states  will  not  be  modeled 
in  the  "suboptimal  filter".  Let  X^  be  the  n^  state  vector  of  the  sub- 
optimal  filter.  The  model  for  X^  is 


( 
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(A-2) 


«k4  + ^ 


Based  upon  the  model  in  (A-2)  (and  possibly  other  constraints) 
a set  of  filter  gains  are  computed.  Consider  the  case  where 
is  a subvector  of  X It  is  of  interest  to  look  at  the  difference 

M a 

between  the  estimate  of  X ^ (denoted  X,  ) and  the  reference  vector 
X^.  In  particular,  consider  the  covariance  matrix  P defined  by 


pk  - ekx,,  - y hx*  - y»Ti 


(A  - 3) 


where  the  * indicates  the  vector  has  been  augmented  with  zeroes 
sufficient  to  make  the  subtraction  well  defined. 

Let  xj  denote  the  vector  difference  (X^  - X^)-  With 
^ j a submatrix  of  ^ ^ and  a submatrix  of  H^,  the  time 
evolution  of  X^  can  be  described  by 


0 


p = <o  p'  <oT  + r rT 

k k, k-1  k-1  k,  k-1  k,  k- 1 k,k-l 


(A  - 4a ) 


‘ W Pkj'ns  - (Kk«k)*jT  + K*\K*T  (A-4b) 


whe  re 


I is  the  n x n identity  matrix 
n s s 


and  the  * indicates  that  the  matrix  with  rows  has  been  aug- 

mented with  (n  - n^)  rows  of  zeroes  so  that  the  matrix 
operations  are  well  defined. 


a 0^0  0^0 

JA  more  general  case  which  includes  errors  in  <p,  T and  H may 
also  be  considered. 
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The  matrices  P,'  and  P,  contain  the  information  about  the 
k k 

variances  of  the  error  in  the  estimate  X^.  With  the  filter  gain 

matrix  Kk  defined,  equations  (A-4a)  and  (A-4b)  are  linear  equations 

for  the  covariance  matrix  P^  for  all  k . It  thus  makes  sense  to 

determine  the  sensitivity  of  the  variance  of  the  filter  states  to  the 

variances  of  the  unmodeled  error  sources.  It  has  been  suggested 

[A-l  ] that  to  do  this,  each  individual  error  source  be  evaluated 

separately  by  propagating  equations  (A-4a)  and  (A-4b)  for  each 

error  source.  However,  for  a certain  class  of  unmodeled  error 

sources,  the  sensitivity  may  be  determined  in  the  presence  of  other 

error  sources  in  just  one  propagation  of  equations  (A-4a)  and  (A-4b). 

0 

To  define  this  class  of  error  sources,  partition  as  follows: 


where 


is  the  subset  of  x£  not  in  the  class  of  interest  Xk  * 

n -vector  where  n * n 
c c i 

is  the  subset  of  which  is  the  class  of  interest 

is  a (n  - n )-vector 

sc 
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where  the  partitioning  is  comformable  with  the  state  vector  parti- 
tion. 

To  show  how  the  sensitivity  of  an  error  source  in  this  class 
may  be  determined  in  the  presence  of  other  error  sources,  the  equa- 
tions for  the  propagation  of  the  inverse  of  will  be  developed.  It 
will  be  apparent  from  these  equations  that  the  effects  of  the  class  of 
error  sources  considered  may  be  easily  removed  from  the  inverse 
covariance  matrix.  Reinverting  the  remaining  part  of  the  inverse 
covariance  matrix  will  yield  the  estimation  error  covariance  matrix 
in  the  absence  of  the  error  source.  The  decrease  in  variance  divided 
by  the  variance  of  the  error  source  used  in  the  covariance  analysis 
defines  the  sensitivity. 

Using  the  state  transition  matrix  and  input  noise  distribution 
matrix  from  (A- 5)  in  (A-4a),  the  following  is  obtained 
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This  may  be  rewritten  as 


where  use  has  been  made  of  the  property  that  ^ j)  = 1 k 

for  any  state  transition  matrix. 

Likewise  the  covariance  after  estimation,  equation  (A-4b) 


may  be  rewritten  as 


!'»  -«wiT  ,A-8» 

s ’ 


I 


181 


Before  proceeding  to  simplify  (A-8)  further,  the  following 
result  which  will  prove  useful  is  noted. 


(A  I' 


1 -a-‘bc-1 


BC  \ 

r1  ) 


(A-9) 


for  arbitrary  size  partitions  when  A and  C are  nonsingular  and  O 
is  a matrix  of  zero. 

Using  the  fact  that  is  the  submatrix  of  H^, which  is  con- 
formable with  the  filter  state  vector  and  (A-9)  in  (A-8), the  following 


is  obtained 


pk  ■ j\  - 


(\  - KA)'1  KRrT  ('„f  • 


■„  -<*W‘ 


(A- 10) 


Equations  (A-7)  and  (A-10)  may  now  be  inverted  with  relative 
ease.  The  only  hard  part  is  the  matrix  sum  sandwiched  in  both 
equations.  Since  they  both  have  the  same  basic  form,  they  can  be 
looked  at  as  inverting  the  following 


M:d)  * Ho) p*' 


(A- 1 1 ) 
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i be  partitioned  as  (using  the  known  symmetry) 
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and  using  the  transpose  of  (A-9) 
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(A- 12) 
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Now  if  in  equation  (A-7)  «>„  ^ k , ^ k , **  = ^ k_,  i. 
identified  with  Q in  equation  (A-12),  then  the  extrapolation  for  the 
inverse  covariance  matrix  is 
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k-1,  k 


Now  in  equation  (A- 10)  identify  with  Q an  n x n matrix 

c c 


(\  - KRK  (\  - KkPk) 


|=  Sk  (A- 14) 


where  the  zeroes  are  matrices  of  appropriate  size  to  make  up  the 

difference  (n  - n,). 

c f 
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Also  let 


"k-K"# 


(A- 15) 


Then  a partition  of  j In 


conformable  with 


the  partition  of 


can  be  obtained  using  equation  (A-9)  as 


Using  equations  (A-10),  (A-12),  (A-14),  and  (A-16),  the  inverse 


of  can  be  expressed 


(A- 17) 
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Now  from  equations  (A- 13)  and  (A- 17)  it  is  apparent  that  the 

upper  left  block  of  the  inverse  covariance  matrix  (i.  e. , the  part 

associated  with  the  state  in  X.  ) is  unaffected  by  the  model  for  the 
2 ~K 

error  states  in  X.  • In  fact,  the  propagation  of  the  upper  left  block 
K 2 

is  the  same  as  it  would  be  if  the  states  X^  were  not  even  considered. 
Thus,  the  effect  of  this  class  of  error  source  may  be  removed  from 
the  covariance  matrix  at  any  point  in  time  simply  by  removing  the 
appropriate  rows  and  columns  from  the  inverse  and  then  reinverting 
the  remaining  submatrix. 

It  may  be  noted  at  this  point  that  the  inverting  of  a large 
covariance  matrix,  elimination  of  a few  rows  and  columns  and  rein- 
version of  the  remaining  matrix  is  a rather  cumbersome  procedure. 
This  is  especially  true  for  large  matrices  and  for  may  sets  of  error 
sources  to  be  examined.  So  at  this  point,  it  will  be  shown  that  the 
desired  result  may  be  obtained  without  inverting  any  large  matrices. 
The  only  inverse  required  has  the  dimension  of  the  error  model  whose 
sensitivity  is  desired. 

From  the  partition  of  P *,  the  desired  matrix  for  the  sensi- 
-1  K 

tivity  analysis  is  Pj  ^ . First  partition  conformable  to  the  parti- 
tion of  P"  * 
k 


(A- 18) 


Then  look  at  the  product  P~*P^ 
blocks  of  (A- 18)  and  in  particular  note 


= I in  terms  of  the  partitioned 
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From  the  relations  in  (A- 19),  the  matrix  P[j  may  be  obtained 
explicitly  in  terms  of  the  submatrices  of  P as  follows.  From  (A- 19a) 


P = (i  _ p \ p"  * 

11  ' P12P12,P11 


(A -20) 


F rom  (A-  19b) 


P = - P P P~ 1 
12  11P12P22 
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With  (A-2 1 ) in  (A-20) 
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So  if  each  error  source  is  considered  individually,  the  only 
inversion  required  is  a matrix  of  the  order  of  the  error  model.  In 
the  case  of  an  error  model  which  contains  only  one  state,  equation 
(A-22)  is  particularly  simple.  In  fact,  if  it  is  only  of  interest  to 
determine  the  variance  reduction,  then  for  each  state  in  X ^ the 
following  holds 


, 2 2 2 

O = <7  (1  - r ) 

i l l 


where 


/ 2 aL  1 

(X.  is  the  new  variance  of  the  i111  state  of  X,  with  the 
l — k 

error  removed 


ctT  is  the  variance  of  the  i^1  state  before  the  error 
i 


is  removed 


r.  is  the  correlation  coefficient  between  the  i state  of 
1 1 

X^  and  the  removed  error  state. 


The  sensitivity  of  the  i^1  state  of  to  the  error  removed 


is  thus 


2 2 

r.O, 

Of.  = -—7- 

1 a2 


where 


a.  is  the  sensitivity  to  the  removed  error  state 


a is  the  variance  of  the  error  state  which  was  used 


to  determine  O.  . 


Reference 

A-l.  Gelb,  A.,  ed.  , Applied  Optimal  Estimation,  the  MIT  Press, 
Cambridge,  Massachusetts,  1974. 
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Appendix  B 

Model  for  Correlated  Ionospheric  Delay  Errors 

The  model  given  in  [1]  for  the  ionospheric  delay  error  indicates 
that  in  addition  to  the  time  correlation,  there  is  a spatial  correlation. 

In  the  case  of  a single  satellite  signal,  this  spatial  correlation  was  con- 
verted to  an  equivalent  time  correlation  by  using  the  constant  speed  of 
the  ionospheric  pierce  point  of  the  signal.  However,  this  spatial  corre- 
lation term  introduces  cross-correlations  among  the  various  satellite 
signals.  Ignoring  these  correlations  in  the  reference  system  model 
would  yield  error  analysis  results  which  are  pessimistically  large. 

A model  which  produces  the  appropriate  autocorrelation  lor  the 
ionospheric  delay  error  is  given  by 


r 


X + 


<Tu(t) 


(B-l) 


where 


CT  is  the  variance  of  the  state  x 
T is  the  associated  time  constant 
u(t)  is  a white  noise  with  unity  PSD 

This  standard  model  for  a Gauss-Markov  process  can  be  used 
to  generate  the  ionospheric  residual  error  for  each  satellite  signal.  In 
addition  to  requirements  on  the  autocorrelation,  it  is  desired  to  have  the 
model  generate  the  proper  cross  correlation.  In  particular,  let  Xj  and 
x^  denote  the  ionospheric  delay  error  states  for  two  satellite  signals. 

*1  ‘ ' T*1  +Vl  Vl'*' 
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From  [1],  the  desired  cross  covariance  is 


R (0)  = <T  , (T,  exp(- tp  /2500  km) 

xx  12 

12 


(B  - 3) 


where 


a.  = « csc[  Ve2  + (18°)2  ] 
i ’ l 


C is  the  correction  residual 


E is  the  elevation  angle  of  the  LOS  vector 
i 


Ap  is  the  distance  between  the  ionospheric  pierce  points 

= |Rj  - R^l  where  R^  is  the  position  of  the  pierce  point 


Using  the  steady  state  values  for  x^  and  x^ 


R , (0)  = E[x  (t)  x,(t)] 
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t t _2t  _€  _U 

= ^d€  ^drj  E[u1(«)u2(17  1)0^2  e Te  Te  T 


Now  let  Uj  and  have  a cross  correlation  of  yfl(t)  , 


2t  t 2« 

- Vze  r * $ • ' T 


de 


°,o2  (f)y 


(B.4) 


Using  (A- 3)  and  (A-4) 


y = (“)  exp(-Ap/2500  km) 

The  required  G matrix  to  generate  the  correlated  driving  noises 
u j and  from  uncorrelated  white  noises  u^  and  u2  can  be  found  from 


[(>*  "2’] 1 e[cQ)("*v  °t] 
s ° -LKv] ci 


• GG 


Then 


GG 


•!* 


ffic2(r)  eitp  <-AP/2500> 


al<r2(r)**p(-Ap/2500)  °\fy 


Th*  extenflion  to  four  correlated  ionospheric  delay  errors  is  now  straight- 
forward. 
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APPENDIX  C 
Glossary  of 

Abbreviations  and  Acronyms 


1 

CANOMIS: 

Covariance  ANalysis  Of  Multisensor  Integrated  System  - 
an  ORINCON  computer  program  for  covariance  analysis 
and  Monte  Carlo  simulation. 

1 

ECI: 

Earth  Centered  Inertial  - a coordinate  frame  with 
origin  at  the  center  of  the  Earth  and  fixed  in  inertial 
space. 

EM-log: 

Electro  Magnetic  log  - a ship's  instrument  to  determine 
speed  with  respect  to  the  water. 

1 

GDE: 

General  Dynamics  Electronics  - the  prime  contractor 
for  the  NAVSTAR  system. 

1 

GDOP: 

Geometric  Dilution  Of  Precision  - a relative  measure 
of  satellite  constellations. 

GPS: 

Global  Positioning  System 

1 

NR: 

Newton  Raphson  - a numerical  technique  for  iterative* 
solution  of  nonlinear  equations. 

PSD: 

Power  Spectral  Density 

t 

RMS: 

Root  Mean  Square  - the  square  root  of  the  average  of 
squares  of  a set  of  numbers. 

RSS: 

Root  Sum  Square  - the  square  root  of  the  sum  of  squares 
of  a set  of  numbers. 

TOA: 

Time-of -Arrival 

WGS: 

World  Geodetic  System 
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